source: palm/trunk/SOURCE/check_parameters.f90 @ 2349

Last change on this file since 2349 was 2348, checked in by kanani, 7 years ago

Parameter check PA0347 added

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/forwind/SOURCE/check_parameters.f901564-1913
File size: 164.8 KB
Line 
1!> @file check_parameters.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
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: check_parameters.f90 2348 2017-08-10 10:40:10Z Giersch $
27! New: Check for simultaneous use of geostrophic wind and u_profile/v_profile
28!
29! 2345 2017-08-09 11:50:30Z Giersch
30! Remove error message PA0156 and the conserve_volume_flow_mode option
31! inflow_profile
32!
33! 2339 2017-08-07 13:55:26Z raasch
34! corrected timestamp in header
35!
36! 2338 2017-08-07 12:15:38Z gronemeier
37! Modularize 1D model
38!
39! 2329 2017-08-03 14:24:56Z knoop
40! Bugfix: index corrected for rho_air and rho_air_zw output
41!
42! 2320 2017-07-21 12:47:43Z suehring
43! Modularize large-scale forcing and nudging
44!
45! 2312 2017-07-14 20:26:51Z hoffmann
46! PA0349 and PA0420 removed.
47!
48! 2300 2017-06-29 13:31:14Z raasch
49! host-specific settings and checks removed
50!
51! 2292 2017-06-20 09:51:42Z schwenkel
52! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
53! includes two more prognostic equations for cloud drop concentration (nc)
54! and cloud water content (qc).
55!
56! 2274 2017-06-09 13:27:48Z Giersch
57! Changed error messages
58!
59! 2271 2017-06-09 12:34:55Z sward
60! roughness-length check altered
61! Error messages fixed
62!
63! 2270 2017-06-09 12:18:47Z maronga
64! Revised numbering (removed 2 timeseries)
65!
66! 2259 2017-06-08 09:09:11Z gronemeier
67! Implemented synthetic turbulence generator
68!
69! 2251 2017-06-06 15:10:46Z Giersch
70!
71! 2250 2017-06-06 15:07:12Z Giersch
72! Doxygen comment added
73!
74! 2248 2017-06-06 13:52:54Z sward
75! Error message fixed
76!
77! 2209 2017-04-19 09:34:46Z kanani
78! Check for plant canopy model output
79!
80! 2200 2017-04-11 11:37:51Z suehring
81! monotonic_adjustment removed
82!
83! 2178 2017-03-17 11:07:39Z hellstea
84! Index limits for perturbations are now set also in case of nested
85! boundary conditions
86!
87! 2169 2017-03-06 18:16:35Z suehring
88! Bugfix, move setting for topography grid convention to init_grid
89!
90! 2118 2017-01-17 16:38:49Z raasch
91! OpenACC related parts of code removed
92!
93! 2088 2016-12-19 16:30:25Z suehring
94! Bugfix in initial salinity profile
95!
96! 2084 2016-12-09 15:59:42Z knoop
97! Removed anelastic + multigrid error message
98!
99! 2052 2016-11-08 15:14:59Z gronemeier
100! Bugfix: remove setting of default value for recycling_width
101!
102! 2050 2016-11-08 15:00:55Z gronemeier
103! Implement turbulent outflow condition
104!
105! 2044 2016-11-02 16:44:25Z knoop
106! Added error code for anelastic approximation
107!
108! 2042 2016-11-02 13:47:31Z suehring
109! Additional checks for wall_heatflux, wall_humidityflux and wall_scalarflux.
110! Bugfix, check for constant_scalarflux.
111!
112! 2040 2016-10-26 16:58:09Z gronemeier
113! Removed check for statistic_regions > 9.
114!
115! 2037 2016-10-26 11:15:40Z knoop
116! Anelastic approximation implemented
117!
118! 2031 2016-10-21 15:11:58Z knoop
119! renamed variable rho to rho_ocean
120!
121! 2026 2016-10-18 10:27:02Z suehring
122! Bugfix, enable output of s*2
123!
124! 2024 2016-10-12 16:42:37Z kanani
125! Added missing CASE, error message and unit for ssws*,
126! increased number of possible output quantities to 500.
127!
128! 2011 2016-09-19 17:29:57Z kanani
129! Flag urban_surface is now defined in module control_parameters,
130! changed prefix for urban surface model output to "usm_",
131! added flag lsf_exception (inipar-Namelist parameter) to allow explicit
132! enabling of large scale forcing together with buildings on flat terrain,
133! introduced control parameter varnamelength for LEN of var.
134!
135! 2007 2016-08-24 15:47:17Z kanani
136! Added checks for the urban surface model,
137! increased counter in DO WHILE loop over data_output (for urban surface output)
138!
139! 2000 2016-08-20 18:09:15Z knoop
140! Forced header and separation lines into 80 columns
141!
142! 1994 2016-08-15 09:52:21Z suehring
143! Add missing check for cloud_physics and cloud_droplets
144!
145! 1992 2016-08-12 15:14:59Z suehring
146! New checks for top_scalarflux
147! Bugfixes concerning data output of profiles in case of passive_scalar
148!
149! 1984 2016-08-01 14:48:05Z suehring
150! Bugfix: checking for bottom and top boundary condition for humidity and scalars
151!
152! 1972 2016-07-26 07:52:02Z maronga
153! Removed check of lai* (done in land surface module)
154!
155! 1970 2016-07-18 14:27:12Z suehring
156! Bugfix, check for ibc_q_b and constant_waterflux in case of land-surface scheme
157!
158! 1962 2016-07-13 09:16:44Z suehring
159! typo removed
160!
161! 1960 2016-07-12 16:34:24Z suehring
162! Separate checks and calculations for humidity and passive scalar. Therefore,
163! a subroutine to calculate vertical gradients is introduced, in order  to reduce
164! redundance.
165! Additional check - large-scale forcing is not implemented for passive scalars
166! so far.
167!
168! 1931 2016-06-10 12:06:59Z suehring
169! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
170!
171! 1929 2016-06-09 16:25:25Z suehring
172! Bugfix in check for use_upstream_for_tke
173!
174! 1918 2016-05-27 14:35:57Z raasch
175! setting of a fixed reference state ('single_value') for ocean runs removed
176!
177! 1914 2016-05-26 14:44:07Z witha
178! Added checks for the wind turbine model
179!
180! 1841 2016-04-07 19:14:06Z raasch
181! redundant location message removed
182!
183! 1833 2016-04-07 14:23:03Z raasch
184! check of spectra quantities moved to spectra_mod
185!
186! 1829 2016-04-07 12:16:29Z maronga
187! Bugfix: output of user defined data required reset of the string 'unit'
188!
189! 1826 2016-04-07 12:01:39Z maronga
190! Moved checks for radiation model to the respective module.
191! Moved checks for plant canopy model to the respective module.
192! Bugfix for check of too large roughness_length
193!
194! 1824 2016-04-07 09:55:29Z gronemeier
195! Check if roughness_length < dz/2. Moved location_message(finished) to the end.
196!
197! 1822 2016-04-07 07:49:42Z hoffmann
198! PALM collision kernel deleted. Collision algorithms are checked for correct
199! spelling.
200!
201! Tails removed.
202! !
203! Checks for use_sgs_for_particles adopted for the use of droplets with
204! use_sgs_for_particles adopted.
205!
206! Unused variables removed.
207!
208! 1817 2016-04-06 15:44:20Z maronga
209! Moved checks for land_surface model to the respective module
210!
211! 1806 2016-04-05 18:55:35Z gronemeier
212! Check for recycling_yshift
213!
214! 1804 2016-04-05 16:30:18Z maronga
215! Removed code for parameter file check (__check)
216!
217! 1795 2016-03-18 15:00:53Z raasch
218! Bugfix: setting of initial scalar profile in ocean
219!
220! 1788 2016-03-10 11:01:04Z maronga
221! Added check for use of most_method = 'lookup' in combination with water
222! surface presribed in the land surface model. Added output of z0q.
223! Syntax layout improved.
224!
225! 1786 2016-03-08 05:49:27Z raasch
226! cpp-direktives for spectra removed, check of spectra level removed
227!
228! 1783 2016-03-06 18:36:17Z raasch
229! netcdf variables and module name changed,
230! check of netcdf precision removed (is done in the netcdf module)
231!
232! 1764 2016-02-28 12:45:19Z raasch
233! output of nest id in run description header,
234! bugfix: check of validity of lateral boundary conditions moved to parin
235!
236! 1762 2016-02-25 12:31:13Z hellstea
237! Introduction of nested domain feature
238!
239! 1745 2016-02-05 13:06:51Z gronemeier
240! Bugfix: check data output intervals to be /= 0.0 in case of parallel NetCDF4
241!
242! 1701 2015-11-02 07:43:04Z maronga
243! Bugfix: definition of rad_net timeseries was missing
244!
245! 1691 2015-10-26 16:17:44Z maronga
246! Added output of Obukhov length (ol) and radiative heating rates for RRTMG.
247! Added checks for use of radiation / lsm with topography.
248!
249! 1682 2015-10-07 23:56:08Z knoop
250! Code annotations made doxygen readable
251!
252! 1606 2015-06-29 10:43:37Z maronga
253! Added check for use of RRTMG without netCDF.
254!
255! 1587 2015-05-04 14:19:01Z maronga
256! Added check for set of albedos when using RRTMG
257!
258! 1585 2015-04-30 07:05:52Z maronga
259! Added support for RRTMG
260!
261! 1575 2015-03-27 09:56:27Z raasch
262! multigrid_fast added as allowed pressure solver
263!
264! 1573 2015-03-23 17:53:37Z suehring
265! Bugfix: check for advection schemes in case of non-cyclic boundary conditions
266!
267! 1557 2015-03-05 16:43:04Z suehring
268! Added checks for monotonic limiter
269!
270! 1555 2015-03-04 17:44:27Z maronga
271! Added output of r_a and r_s. Renumbering of LSM PA-messages.
272!
273! 1553 2015-03-03 17:33:54Z maronga
274! Removed check for missing soil temperature profile as default values were added.
275!
276! 1551 2015-03-03 14:18:16Z maronga
277! Added various checks for land surface and radiation model. In the course of this
278! action, the length of the variable var has to be increased
279!
280! 1504 2014-12-04 09:23:49Z maronga
281! Bugfix: check for large_scale forcing got lost.
282!
283! 1500 2014-12-03 17:42:41Z maronga
284! Boundary conditions changed to dirichlet for land surface model
285!
286! 1496 2014-12-02 17:25:50Z maronga
287! Added checks for the land surface model
288!
289! 1484 2014-10-21 10:53:05Z kanani
290! Changes due to new module structure of the plant canopy model:
291!   module plant_canopy_model_mod added,
292!   checks regarding usage of new method for leaf area density profile
293!   construction added,
294!   lad-profile construction moved to new subroutine init_plant_canopy within
295!   the module plant_canopy_model_mod,
296!   drag_coefficient renamed to canopy_drag_coeff.
297! Missing KIND-attribute for REAL constant added
298!
299! 1455 2014-08-29 10:47:47Z heinze
300! empty time records in volume, cross-section and masked data output prevented
301! in case of non-parallel netcdf-output in restart runs
302!
303! 1429 2014-07-15 12:53:45Z knoop
304! run_description_header exended to provide ensemble_member_nr if specified
305!
306! 1425 2014-07-05 10:57:53Z knoop
307! bugfix: perturbation domain modified for parallel random number generator
308!
309! 1402 2014-05-09 14:25:13Z raasch
310! location messages modified
311!
312! 1400 2014-05-09 14:03:54Z knoop
313! Check random generator extended by option random-parallel
314!
315! 1384 2014-05-02 14:31:06Z raasch
316! location messages added
317!
318! 1365 2014-04-22 15:03:56Z boeske
319! Usage of large scale forcing for pt and q enabled:
320! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
321! large scale subsidence (td_sub_lpt, td_sub_q)
322! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
323! check if use_subsidence_tendencies is used correctly
324!
325! 1361 2014-04-16 15:17:48Z hoffmann
326! PA0363 removed
327! PA0362 changed
328!
329! 1359 2014-04-11 17:15:14Z hoffmann
330! Do not allow the execution of PALM with use_particle_tails, since particle
331! tails are currently not supported by our new particle structure.
332!
333! PA0084 not necessary for new particle structure
334!
335! 1353 2014-04-08 15:21:23Z heinze
336! REAL constants provided with KIND-attribute
337!
338! 1330 2014-03-24 17:29:32Z suehring
339! In case of SGS-particle velocity advection of TKE is also allowed with
340! dissipative 5th-order scheme.
341!
342! 1327 2014-03-21 11:00:16Z raasch
343! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
344! bugfix: duplicate error message 56 removed,
345! check of data_output_format and do3d_compress removed
346!
347! 1322 2014-03-20 16:38:49Z raasch
348! some REAL constants defined as wp-kind
349!
350! 1320 2014-03-20 08:40:49Z raasch
351! Kind-parameters added to all INTEGER and REAL declaration statements,
352! kinds are defined in new module kinds,
353! revision history before 2012 removed,
354! comment fields (!:) to be used for variable explanations added to
355! all variable declaration statements
356!
357! 1308 2014-03-13 14:58:42Z fricke
358! +netcdf_data_format_save
359! Calculate fixed number of output time levels for parallel netcdf output.
360! For masked data, parallel netcdf output is not tested so far, hence
361! netcdf_data_format is switched back to non-paralell output.
362!
363! 1299 2014-03-06 13:15:21Z heinze
364! enable usage of large_scale subsidence in combination with large_scale_forcing
365! output for profile of large scale vertical velocity w_subs added
366!
367! 1276 2014-01-15 13:40:41Z heinze
368! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
369!
370! 1241 2013-10-30 11:36:58Z heinze
371! output for profiles of ug and vg added
372! set surface_heatflux and surface_waterflux also in dependence on
373! large_scale_forcing
374! checks for nudging and large scale forcing from external file
375!
376! 1236 2013-09-27 07:21:13Z raasch
377! check number of spectra levels
378!
379! 1216 2013-08-26 09:31:42Z raasch
380! check for transpose_compute_overlap (temporary)
381!
382! 1214 2013-08-21 12:29:17Z kanani
383! additional check for simultaneous use of vertical grid stretching
384! and particle advection
385!
386! 1212 2013-08-15 08:46:27Z raasch
387! checks for poisfft_hybrid removed
388!
389! 1210 2013-08-14 10:58:20Z raasch
390! check for fftw
391!
392! 1179 2013-06-14 05:57:58Z raasch
393! checks and settings of buoyancy parameters and switches revised,
394! initial profile for rho_ocean added to hom (id=77)
395!
396! 1174 2013-05-31 10:28:08Z gryschka
397! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
398!
399! 1159 2013-05-21 11:58:22Z fricke
400! bc_lr/ns_dirneu/neudir removed
401!
402! 1115 2013-03-26 18:16:16Z hoffmann
403! unused variables removed
404! drizzle can be used without precipitation
405!
406! 1111 2013-03-08 23:54:10Z raasch
407! ibc_p_b = 2 removed
408!
409! 1103 2013-02-20 02:15:53Z raasch
410! Bugfix: turbulent inflow must not require cyclic fill in restart runs
411!
412! 1092 2013-02-02 11:24:22Z raasch
413! unused variables removed
414!
415! 1069 2012-11-28 16:18:43Z maronga
416! allow usage of topography in combination with cloud physics
417!
418! 1065 2012-11-22 17:42:36Z hoffmann
419! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
420!         precipitation in order to save computational resources.
421!
422! 1060 2012-11-21 07:19:51Z raasch
423! additional check for parameter turbulent_inflow
424!
425! 1053 2012-11-13 17:11:03Z hoffmann
426! necessary changes for the new two-moment cloud physics scheme added:
427! - check cloud physics scheme (Kessler or Seifert and Beheng)
428! - plant_canopy is not allowed
429! - currently, only cache loop_optimization is allowed
430! - initial profiles of nr, qr
431! - boundary condition of nr, qr
432! - check output quantities (qr, nr, prr)
433!
434! 1036 2012-10-22 13:43:42Z raasch
435! code put under GPL (PALM 3.9)
436!
437! 1031/1034 2012-10-22 11:32:49Z raasch
438! check of netcdf4 parallel file support
439!
440! 1019 2012-09-28 06:46:45Z raasch
441! non-optimized version of prognostic_equations not allowed any more
442!
443! 1015 2012-09-27 09:23:24Z raasch
444! acc allowed for loop optimization,
445! checks for adjustment of mixing length to the Prandtl mixing length removed
446!
447! 1003 2012-09-14 14:35:53Z raasch
448! checks for cases with unequal subdomain sizes removed
449!
450! 1001 2012-09-13 14:08:46Z raasch
451! all actions concerning leapfrog- and upstream-spline-scheme removed
452!
453! 996 2012-09-07 10:41:47Z raasch
454! little reformatting
455
456! 978 2012-08-09 08:28:32Z fricke
457! setting of bc_lr/ns_dirneu/neudir
458! outflow damping layer removed
459! check for z0h*
460! check for pt_damping_width
461!
462! 964 2012-07-26 09:14:24Z raasch
463! check of old profil-parameters removed
464!
465! 940 2012-07-09 14:31:00Z raasch
466! checks for parameter neutral
467!
468! 924 2012-06-06 07:44:41Z maronga
469! Bugfix: preprocessor directives caused error during compilation
470!
471! 892 2012-05-02 13:51:44Z maronga
472! Bugfix for parameter file check ( excluding __netcdf4 )
473!
474! 866 2012-03-28 06:44:41Z raasch
475! use only 60% of the geostrophic wind as translation speed in case of Galilean
476! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
477! timestep
478!
479! 861 2012-03-26 14:18:34Z suehring
480! Check for topography and ws-scheme removed.
481! Check for loop_optimization = 'vector' and ws-scheme removed.
482!
483! 845 2012-03-07 10:23:05Z maronga
484! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
485!
486! 828 2012-02-21 12:00:36Z raasch
487! check of collision_kernel extended
488!
489! 825 2012-02-19 03:03:44Z raasch
490! check for collision_kernel and curvature_solution_effects
491!
492! 809 2012-01-30 13:32:58Z maronga
493! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
494!
495! 807 2012-01-25 11:53:51Z maronga
496! New cpp directive "__check" implemented which is used by check_namelist_files
497!
498! Revision 1.1  1997/08/26 06:29:23  raasch
499! Initial revision
500!
501!
502! Description:
503! ------------
504!> Check control parameters and deduce further quantities.
505!------------------------------------------------------------------------------!
506 SUBROUTINE check_parameters
507
508
509    USE arrays_3d
510    USE cloud_parameters
511    USE constants
512    USE control_parameters
513    USE dvrp_variables
514    USE grid_variables
515    USE indices
516    USE land_surface_model_mod,                                                &
517        ONLY: lsm_check_data_output, lsm_check_data_output_pr,                 &
518              lsm_check_parameters
519
520    USE lsf_nudging_mod,                                                       &
521        ONLY:  lsf_nudging_check_parameters, lsf_nudging_check_data_output_pr, &
522               qsws_surf, shf_surf
523
524    USE kinds
525
526    USE model_1d_mod,                                                          &
527        ONLY:  damp_level_1d, damp_level_ind_1d
528
529    USE netcdf_interface,                                                      &
530        ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
531               netcdf_data_format_string, dots_unit, heatflux_output_unit,     &
532               waterflux_output_unit, momentumflux_output_unit
533
534    USE particle_attributes
535    USE pegrid
536    USE plant_canopy_model_mod,                                                &
537        ONLY:  pcm_check_data_output, pcm_check_parameters, plant_canopy
538
539    USE pmc_interface,                                                         &
540        ONLY:  cpl_id, nested_run
541
542    USE profil_parameter
543    USE radiation_model_mod,                                                   &
544        ONLY:  radiation, radiation_check_data_output,                         &
545               radiation_check_data_output_pr, radiation_check_parameters
546    USE spectra_mod,                                                           &
547        ONLY:  calculate_spectra, spectra_check_parameters
548    USE statistics
549    USE subsidence_mod
550    USE statistics
551    USE synthetic_turbulence_generator_mod,                                    &
552        ONLY:  stg_check_parameters
553    USE transpose_indices
554    USE urban_surface_mod,                                                     &
555        ONLY:  usm_check_data_output, usm_check_parameters
556    USE wind_turbine_model_mod,                                                &
557        ONLY:  wtm_check_parameters, wind_turbine
558
559
560    IMPLICIT NONE
561
562    CHARACTER (LEN=1)   ::  sq                       !<
563    CHARACTER (LEN=varnamelength)  ::  var           !<
564    CHARACTER (LEN=7)   ::  unit                     !<
565    CHARACTER (LEN=8)   ::  date                     !<
566    CHARACTER (LEN=10)  ::  time                     !<
567    CHARACTER (LEN=20)  ::  ensemble_string          !<
568    CHARACTER (LEN=15)  ::  nest_string              !<
569    CHARACTER (LEN=40)  ::  coupling_string          !<
570    CHARACTER (LEN=100) ::  action                   !<
571
572    INTEGER(iwp) ::  i                               !<
573    INTEGER(iwp) ::  ilen                            !<
574    INTEGER(iwp) ::  j                               !<
575    INTEGER(iwp) ::  k                               !<
576    INTEGER(iwp) ::  kk                              !<
577    INTEGER(iwp) ::  netcdf_data_format_save         !<
578    INTEGER(iwp) ::  position                        !<
579
580    LOGICAL     ::  found                            !<
581
582    REAL(wp)    ::  dum                              !<
583    REAL(wp)    ::  gradient                         !<
584    REAL(wp)    ::  remote = 0.0_wp                  !<
585    REAL(wp)    ::  simulation_time_since_reference  !<
586
587
588    CALL location_message( 'checking parameters', .FALSE. )
589
590!
591!-- Check for overlap combinations, which are not realized yet
592    IF ( transpose_compute_overlap )  THEN
593       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
594    ENDIF
595
596!
597!-- Check the coupling mode
598!> @todo Check if any queries for other coupling modes (e.g. precursor_ocean) are missing
599    IF ( coupling_mode /= 'uncoupled'            .AND.  &
600         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
601         coupling_mode /= 'ocean_to_atmosphere' )  THEN
602       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
603       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
604    ENDIF
605
606!
607!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
608    IF ( coupling_mode /= 'uncoupled')  THEN
609
610       IF ( dt_coupling == 9999999.9_wp )  THEN
611          message_string = 'dt_coupling is not set but required for coup' //   &
612                           'ling mode "' //  TRIM( coupling_mode ) // '"'
613          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
614       ENDIF
615
616#if defined( __parallel )
617
618!
619!--    NOTE: coupled runs have not been implemented in the check_namelist_files
620!--    program.
621!--    check_namelist_files will need the following information of the other
622!--    model (atmosphere/ocean).
623!       dt_coupling = remote
624!       dt_max = remote
625!       restart_time = remote
626!       dt_restart= remote
627!       simulation_time_since_reference = remote
628!       dx = remote
629
630
631       IF ( myid == 0 ) THEN
632          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter,  &
633                         ierr )
634          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
635                         status, ierr )
636       ENDIF
637       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
638
639       IF ( dt_coupling /= remote )  THEN
640          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
641                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
642                 'dt_coupling_remote = ', remote
643          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
644       ENDIF
645       IF ( dt_coupling <= 0.0_wp )  THEN
646
647          IF ( myid == 0  ) THEN
648             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
649             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
650                            status, ierr )
651          ENDIF
652          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
653
654          dt_coupling = MAX( dt_max, remote )
655          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
656                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
657                 'MAX(dt_max(A,O)) = ', dt_coupling
658          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
659       ENDIF
660
661       IF ( myid == 0 ) THEN
662          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
663                         ierr )
664          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
665                         status, ierr )
666       ENDIF
667       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
668
669       IF ( restart_time /= remote )  THEN
670          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
671                 '": restart_time = ', restart_time, '& is not equal to ',     &
672                 'restart_time_remote = ', remote
673          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
674       ENDIF
675
676       IF ( myid == 0 ) THEN
677          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter,   &
678                         ierr )
679          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
680                         status, ierr )
681       ENDIF
682       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
683
684       IF ( dt_restart /= remote )  THEN
685          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
686                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
687                 'dt_restart_remote = ', remote
688          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
689       ENDIF
690
691       simulation_time_since_reference = end_time - coupling_start_time
692
693       IF  ( myid == 0 ) THEN
694          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
695                         14, comm_inter, ierr )
696          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
697                         status, ierr )
698       ENDIF
699       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
700
701       IF ( simulation_time_since_reference /= remote )  THEN
702          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
703                 '": simulation_time_since_reference = ',                      &
704                 simulation_time_since_reference, '& is not equal to ',        &
705                 'simulation_time_since_reference_remote = ', remote
706          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
707       ENDIF
708
709       IF ( myid == 0 ) THEN
710          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
711          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
712                                                             status, ierr )
713       ENDIF
714       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
715
716
717       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
718
719          IF ( dx < remote ) THEN
720             WRITE( message_string, * ) 'coupling mode "',                     &
721                   TRIM( coupling_mode ),                                      &
722           '": dx in Atmosphere is not equal to or not larger than dx in ocean'
723             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
724          ENDIF
725
726          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
727             WRITE( message_string, * ) 'coupling mode "',                     &
728                    TRIM( coupling_mode ),                                     &
729             '": Domain size in x-direction is not equal in ocean and atmosphere'
730             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
731          ENDIF
732
733       ENDIF
734
735       IF ( myid == 0) THEN
736          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
737          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter,       &
738                         status, ierr )
739       ENDIF
740       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
741
742       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
743
744          IF ( dy < remote )  THEN
745             WRITE( message_string, * ) 'coupling mode "',                     &
746                    TRIM( coupling_mode ),                                     &
747                 '": dy in Atmosphere is not equal to or not larger than dy in ocean'
748             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
749          ENDIF
750
751          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
752             WRITE( message_string, * ) 'coupling mode "',                     &
753                   TRIM( coupling_mode ),                                      &
754             '": Domain size in y-direction is not equal in ocean and atmosphere'
755             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
756          ENDIF
757
758          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
759             WRITE( message_string, * ) 'coupling mode "',                     &
760                   TRIM( coupling_mode ),                                      &
761             '": nx+1 in ocean is not divisible by nx+1 in', &
762             ' atmosphere without remainder'
763             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
764          ENDIF
765
766          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
767             WRITE( message_string, * ) 'coupling mode "',                     &
768                   TRIM( coupling_mode ),                                      &
769             '": ny+1 in ocean is not divisible by ny+1 in', &
770             ' atmosphere without remainder'
771
772             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
773          ENDIF
774
775       ENDIF
776#else
777       WRITE( message_string, * ) 'coupling requires PALM to be called with',  &
778            ' ''mrun -K parallel'''
779       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
780#endif
781    ENDIF
782
783#if defined( __parallel )
784!
785!-- Exchange via intercommunicator
786    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
787       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter,     &
788                      ierr )
789    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
790       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19,          &
791                      comm_inter, status, ierr )
792    ENDIF
793    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
794
795#endif
796
797
798!
799!-- Generate the file header which is used as a header for most of PALM's
800!-- output files
801    CALL DATE_AND_TIME( date, time )
802    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
803    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
804    IF ( coupling_mode == 'uncoupled' )  THEN
805       coupling_string = ''
806    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
807       coupling_string = ' coupled (atmosphere)'
808    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
809       coupling_string = ' coupled (ocean)'
810    ENDIF
811    IF ( ensemble_member_nr /= 0 )  THEN
812       WRITE( ensemble_string, '(2X,A,I2.2)' )  'en-no: ', ensemble_member_nr
813    ELSE
814       ensemble_string = ''
815    ENDIF
816    IF ( nested_run )  THEN
817       WRITE( nest_string, '(2X,A,I2.2)' )  'nest-id: ', cpl_id
818    ELSE
819       nest_string = ''
820    ENDIF
821
822    WRITE ( run_description_header,                                            &
823            '(A,2X,A,2X,A,A,A,I2.2,A,A,A,2X,A,A,2X,A,1X,A)' )                  &
824          TRIM( version ), TRIM( revision ), 'run: ',                          &
825          TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),         &
826          TRIM( nest_string ), TRIM( ensemble_string), 'host: ', TRIM( host ), &
827          run_date, run_time
828
829!
830!-- Check the general loop optimization method
831    SELECT CASE ( TRIM( loop_optimization ) )
832
833       CASE ( 'cache', 'vector' )
834          CONTINUE
835
836       CASE DEFAULT
837          message_string = 'illegal value given for loop_optimization: "' //   &
838                           TRIM( loop_optimization ) // '"'
839          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
840
841    END SELECT
842
843!
844!-- Check if vertical grid stretching is used together with particles
845    IF ( dz_stretch_level < 100000.0_wp .AND. particle_advection )  THEN
846       message_string = 'Vertical grid stretching is not allowed together ' // &
847                        'with particle advection.'
848       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
849    ENDIF
850
851!
852!-- Check topography setting (check for illegal parameter combinations)
853    IF ( topography /= 'flat' )  THEN
854       action = ' '
855       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'      &
856          )  THEN
857          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
858       ENDIF
859       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )&
860       THEN
861          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
862       ENDIF
863       IF ( psolver == 'sor' )  THEN
864          WRITE( action, '(A,A)' )  'psolver = ', psolver
865       ENDIF
866       IF ( sloping_surface )  THEN
867          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
868       ENDIF
869       IF ( galilei_transformation )  THEN
870          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
871       ENDIF
872       IF ( cloud_physics )  THEN
873          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
874       ENDIF
875       IF ( cloud_droplets )  THEN
876          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
877       ENDIF
878       IF ( .NOT. constant_flux_layer )  THEN
879          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
880       ENDIF
881       IF ( action /= ' ' )  THEN
882          message_string = 'a non-flat topography does not allow ' //          &
883                           TRIM( action )
884          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
885       ENDIF
886
887    ENDIF
888
889!
890!-- Check approximation
891    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.   &
892         TRIM( approximation ) /= 'anelastic' )  THEN
893       message_string = 'unknown approximation: approximation = "' //    &
894                        TRIM( approximation ) // '"'
895       CALL message( 'check_parameters', 'PA0446', 1, 2, 0, 6, 0 )
896    ENDIF
897
898!
899!-- Check approximation requirements
900    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
901         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
902       message_string = 'Anelastic approximation requires: ' //                &
903                        'momentum_advec = "ws-scheme"'
904       CALL message( 'check_parameters', 'PA0447', 1, 2, 0, 6, 0 )
905    ENDIF
906    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
907         TRIM( psolver ) == 'multigrid' )  THEN
908       message_string = 'Anelastic approximation currently only supports: ' // &
909                        'psolver = "poisfft", ' // &
910                        'psolver = "sor" and ' // &
911                        'psolver = "multigrid_noopt"'
912       CALL message( 'check_parameters', 'PA0448', 1, 2, 0, 6, 0 )
913    ENDIF
914    IF ( TRIM( approximation ) == 'anelastic'   .AND.   &
915         conserve_volume_flow )  THEN
916       message_string = 'Anelastic approximation is not allowed with:' //      &
917                        'conserve_volume_flow = .TRUE.'
918       CALL message( 'check_parameters', 'PA0449', 1, 2, 0, 6, 0 )
919    ENDIF
920
921!
922!-- Check flux input mode
923    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.   &
924         TRIM( flux_input_mode ) /= 'kinematic'  .AND.   &
925         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
926       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
927                        TRIM( flux_input_mode ) // '"'
928       CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
929    ENDIF
930!-- Set flux input mode according to approximation if applicable
931    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
932       IF ( TRIM( approximation ) == 'anelastic' )  THEN
933          flux_input_mode = 'dynamic'
934       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
935          flux_input_mode = 'kinematic'
936       ENDIF
937    ENDIF
938
939!
940!-- Check flux output mode
941    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.   &
942         TRIM( flux_output_mode ) /= 'kinematic'  .AND.   &
943         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
944       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
945                        TRIM( flux_output_mode ) // '"'
946       CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
947    ENDIF
948!-- Set flux output mode according to approximation if applicable
949    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
950       IF ( TRIM( approximation ) == 'anelastic' )  THEN
951          flux_output_mode = 'dynamic'
952       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
953          flux_output_mode = 'kinematic'
954       ENDIF
955    ENDIF
956
957
958!
959!-- When the land surface model is used, the flux output must be dynamic.
960    IF ( land_surface )  THEN
961       flux_output_mode = 'dynamic'
962    ENDIF
963
964!
965!-- set the flux output units according to flux_output_mode
966    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
967        heatflux_output_unit              = 'K m/s'
968        waterflux_output_unit             = 'kg/kg m/s'
969        momentumflux_output_unit          = 'm2/s2'
970    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
971        heatflux_output_unit              = 'W/m2'
972        waterflux_output_unit             = 'W/m2'
973        momentumflux_output_unit          = 'N/m2'
974    ENDIF
975
976!-- set time series output units for fluxes
977    dots_unit(14:16) = heatflux_output_unit
978    dots_unit(21)    = waterflux_output_unit
979    dots_unit(19:20) = momentumflux_output_unit
980
981!
982!-- Check ocean setting
983    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
984         TRIM( coupling_char ) == '_O' .AND.                                   &
985         .NOT. ocean )  THEN
986
987!
988!--    Check whether an (uncoupled) atmospheric run has been declared as an
989!--    ocean run (this setting is done via mrun-option -y)
990       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
991                        TRIM( coupling_char ) // '" set by mrun-option "-y"'
992       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
993
994    ENDIF
995!
996!-- Check cloud scheme
997    IF ( cloud_scheme == 'saturation_adjust' )  THEN
998       microphysics_sat_adjust = .TRUE.
999       microphysics_seifert    = .FALSE.
1000       microphysics_kessler    = .FALSE.
1001       precipitation           = .FALSE.
1002    ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
1003       microphysics_sat_adjust = .FALSE.
1004       microphysics_seifert    = .TRUE.
1005       microphysics_kessler    = .FALSE.
1006       microphysics_morrison  = .FALSE.
1007       precipitation           = .TRUE.
1008    ELSEIF ( cloud_scheme == 'kessler' )  THEN
1009       microphysics_sat_adjust = .FALSE.
1010       microphysics_seifert    = .FALSE.
1011       microphysics_kessler    = .TRUE.
1012       microphysics_morrison   = .FALSE.
1013       precipitation           = .TRUE.
1014    ELSEIF ( cloud_scheme == 'morrison' )  THEN
1015       microphysics_sat_adjust = .FALSE.
1016       microphysics_seifert    = .TRUE.
1017       microphysics_kessler    = .FALSE.
1018       microphysics_morrison   = .TRUE.
1019       precipitation           = .TRUE.
1020    ELSE
1021       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
1022                        TRIM( cloud_scheme ) // '"'
1023       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
1024    ENDIF
1025!
1026!-- Check whether there are any illegal values
1027!-- Pressure solver:
1028    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
1029         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_noopt' )  THEN
1030       message_string = 'unknown solver for perturbation pressure: psolver' // &
1031                        ' = "' // TRIM( psolver ) // '"'
1032       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
1033    ENDIF
1034
1035    IF ( psolver(1:9) == 'multigrid' )  THEN
1036       IF ( cycle_mg == 'w' )  THEN
1037          gamma_mg = 2
1038       ELSEIF ( cycle_mg == 'v' )  THEN
1039          gamma_mg = 1
1040       ELSE
1041          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
1042                           TRIM( cycle_mg ) // '"'
1043          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
1044       ENDIF
1045    ENDIF
1046
1047    IF ( fft_method /= 'singleton-algorithm'  .AND.                            &
1048         fft_method /= 'temperton-algorithm'  .AND.                            &
1049         fft_method /= 'fftw'                 .AND.                            &
1050         fft_method /= 'system-specific' )  THEN
1051       message_string = 'unknown fft-algorithm: fft_method = "' //             &
1052                        TRIM( fft_method ) // '"'
1053       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
1054    ENDIF
1055
1056    IF( momentum_advec == 'ws-scheme' .AND.                                    &
1057        .NOT. call_psolver_at_all_substeps  ) THEN
1058        message_string = 'psolver must be called at each RK3 substep when "'// &
1059                      TRIM(momentum_advec) // ' "is used for momentum_advec'
1060        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
1061    END IF
1062!
1063!-- Advection schemes:
1064    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
1065    THEN
1066       message_string = 'unknown advection scheme: momentum_advec = "' //      &
1067                        TRIM( momentum_advec ) // '"'
1068       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
1069    ENDIF
1070    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
1071           .AND. ( timestep_scheme == 'euler' .OR.                             &
1072                   timestep_scheme == 'runge-kutta-2' ) )                      &
1073    THEN
1074       message_string = 'momentum_advec or scalar_advec = "'                   &
1075         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
1076         TRIM( timestep_scheme ) // '"'
1077       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
1078    ENDIF
1079    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
1080         scalar_advec /= 'bc-scheme' )                                         &
1081    THEN
1082       message_string = 'unknown advection scheme: scalar_advec = "' //        &
1083                        TRIM( scalar_advec ) // '"'
1084       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
1085    ENDIF
1086    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
1087    THEN
1088       message_string = 'advection_scheme scalar_advec = "'                    &
1089         // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
1090         TRIM( loop_optimization ) // '"'
1091       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
1092    ENDIF
1093
1094    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.             &
1095         .NOT. use_upstream_for_tke  .AND.                                     &
1096         scalar_advec /= 'ws-scheme'                                           &
1097       )  THEN
1098       use_upstream_for_tke = .TRUE.
1099       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //    &
1100                        'use_sgs_for_particles = .TRUE. '          //          &
1101                        'and scalar_advec /= ws-scheme'
1102       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
1103    ENDIF
1104
1105!
1106!-- Set LOGICAL switches to enhance performance
1107    IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
1108    IF ( scalar_advec   == 'ws-scheme' )  ws_scheme_sca = .TRUE.
1109
1110
1111!
1112!-- Timestep schemes:
1113    SELECT CASE ( TRIM( timestep_scheme ) )
1114
1115       CASE ( 'euler' )
1116          intermediate_timestep_count_max = 1
1117
1118       CASE ( 'runge-kutta-2' )
1119          intermediate_timestep_count_max = 2
1120
1121       CASE ( 'runge-kutta-3' )
1122          intermediate_timestep_count_max = 3
1123
1124       CASE DEFAULT
1125          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
1126                           TRIM( timestep_scheme ) // '"'
1127          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
1128
1129    END SELECT
1130
1131    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
1132         .AND. timestep_scheme(1:5) == 'runge' ) THEN
1133       message_string = 'momentum advection scheme "' // &
1134                        TRIM( momentum_advec ) // '" & does not work with ' // &
1135                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
1136       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
1137    ENDIF
1138!
1139!-- Check for proper settings for microphysics
1140    IF ( cloud_physics  .AND.  cloud_droplets )  THEN
1141       message_string = 'cloud_physics = .TRUE. is not allowed with ' //  &
1142                        'cloud_droplets = .TRUE.'
1143       CALL message( 'check_parameters', 'PA0442', 1, 2, 0, 6, 0 )
1144    ENDIF
1145
1146!
1147!-- Collision kernels:
1148    SELECT CASE ( TRIM( collision_kernel ) )
1149
1150       CASE ( 'hall', 'hall_fast' )
1151          hall_kernel = .TRUE.
1152
1153       CASE ( 'wang', 'wang_fast' )
1154          wang_kernel = .TRUE.
1155
1156       CASE ( 'none' )
1157
1158
1159       CASE DEFAULT
1160          message_string = 'unknown collision kernel: collision_kernel = "' // &
1161                           TRIM( collision_kernel ) // '"'
1162          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
1163
1164    END SELECT
1165    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
1166
1167    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1168         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1169!
1170!--    No restart run: several initialising actions are possible
1171       action = initializing_actions
1172       DO  WHILE ( TRIM( action ) /= '' )
1173          position = INDEX( action, ' ' )
1174          SELECT CASE ( action(1:position-1) )
1175
1176             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
1177                    'by_user', 'initialize_vortex', 'initialize_ptanom' )
1178                action = action(position+1:)
1179
1180             CASE DEFAULT
1181                message_string = 'initializing_action = "' //                  &
1182                                 TRIM( action ) // '" unkown or not allowed'
1183                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
1184
1185          END SELECT
1186       ENDDO
1187    ENDIF
1188
1189    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
1190         conserve_volume_flow ) THEN
1191         message_string = 'initializing_actions = "initialize_vortex"' //      &
1192                        ' is not allowed with conserve_volume_flow = .T.'
1193       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
1194    ENDIF
1195
1196
1197    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1198         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1199       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1200                        ' and "set_1d-model_profiles" are not allowed ' //     &
1201                        'simultaneously'
1202       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
1203    ENDIF
1204
1205    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1206         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
1207       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1208                        ' and "by_user" are not allowed simultaneously'
1209       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
1210    ENDIF
1211
1212    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                  &
1213         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1214       message_string = 'initializing_actions = "by_user" and ' //             &
1215                        '"set_1d-model_profiles" are not allowed simultaneously'
1216       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
1217    ENDIF
1218
1219    IF ( cloud_physics  .AND.  .NOT.  humidity )  THEN
1220       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ',   &
1221              'not allowed with humidity = ', humidity
1222       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
1223    ENDIF
1224
1225    IF ( humidity  .AND.  sloping_surface )  THEN
1226       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //   &
1227                        'are not allowed simultaneously'
1228       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
1229    ENDIF
1230
1231!
1232!-- When synthetic turbulence generator is used, peform addtional checks
1233    IF ( synthetic_turbulence_generator )  CALL stg_check_parameters
1234!
1235!-- When plant canopy model is used, peform addtional checks
1236    IF ( plant_canopy )  CALL pcm_check_parameters
1237
1238!
1239!-- Additional checks for spectra
1240    IF ( calculate_spectra )  CALL spectra_check_parameters
1241!
1242!-- When land surface model is used, perform additional checks
1243    IF ( land_surface )  CALL lsm_check_parameters
1244
1245!
1246!-- When urban surface model is used, perform additional checks
1247    IF ( urban_surface )  CALL usm_check_parameters
1248
1249!
1250!-- If wind turbine model is used, perform additional checks
1251    IF ( wind_turbine )  CALL wtm_check_parameters
1252!
1253!-- When radiation model is used, peform addtional checks
1254    IF ( radiation )  CALL radiation_check_parameters
1255!
1256!-- When large-scale forcing or nudging is used, peform addtional checks
1257    IF ( large_scale_forcing  .OR.  nudging )  CALL lsf_nudging_check_parameters
1258
1259
1260
1261    IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
1262                 loop_optimization == 'vector' )                               &
1263         .AND.  cloud_physics  .AND.  microphysics_seifert )  THEN
1264       message_string = 'cloud_scheme = seifert_beheng requires ' //           &
1265                        'loop_optimization = "cache" or "vector"'
1266       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
1267    ENDIF
1268
1269!
1270!-- In case of no model continuation run, check initialising parameters and
1271!-- deduce further quantities
1272    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1273
1274!
1275!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
1276       pt_init = pt_surface
1277       IF ( humidity       )  q_init  = q_surface
1278       IF ( ocean          )  sa_init = sa_surface
1279       IF ( passive_scalar )  s_init  = s_surface
1280
1281!
1282!--
1283!--    If required, compute initial profile of the geostrophic wind
1284!--    (component ug)
1285       i = 1
1286       gradient = 0.0_wp
1287
1288       IF (  .NOT.  ocean )  THEN
1289
1290          ug_vertical_gradient_level_ind(1) = 0
1291          ug(0) = ug_surface
1292          DO  k = 1, nzt+1
1293             IF ( i < 11 )  THEN
1294                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.              &
1295                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
1296                   gradient = ug_vertical_gradient(i) / 100.0_wp
1297                   ug_vertical_gradient_level_ind(i) = k - 1
1298                   i = i + 1
1299                ENDIF
1300             ENDIF
1301             IF ( gradient /= 0.0_wp )  THEN
1302                IF ( k /= 1 )  THEN
1303                   ug(k) = ug(k-1) + dzu(k) * gradient
1304                ELSE
1305                   ug(k) = ug_surface + dzu(k) * gradient
1306                ENDIF
1307             ELSE
1308                ug(k) = ug(k-1)
1309             ENDIF
1310          ENDDO
1311
1312       ELSE
1313
1314          ug_vertical_gradient_level_ind(1) = nzt+1
1315          ug(nzt+1) = ug_surface
1316          DO  k = nzt, nzb, -1
1317             IF ( i < 11 )  THEN
1318                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.              &
1319                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
1320                   gradient = ug_vertical_gradient(i) / 100.0_wp
1321                   ug_vertical_gradient_level_ind(i) = k + 1
1322                   i = i + 1
1323                ENDIF
1324             ENDIF
1325             IF ( gradient /= 0.0_wp )  THEN
1326                IF ( k /= nzt )  THEN
1327                   ug(k) = ug(k+1) - dzu(k+1) * gradient
1328                ELSE
1329                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
1330                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
1331                ENDIF
1332             ELSE
1333                ug(k) = ug(k+1)
1334             ENDIF
1335          ENDDO
1336
1337       ENDIF
1338
1339!
1340!--    In case of no given gradients for ug, choose a zero gradient
1341       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1342          ug_vertical_gradient_level(1) = 0.0_wp
1343       ENDIF
1344
1345!
1346!--
1347!--    If required, compute initial profile of the geostrophic wind
1348!--    (component vg)
1349       i = 1
1350       gradient = 0.0_wp
1351
1352       IF (  .NOT.  ocean )  THEN
1353
1354          vg_vertical_gradient_level_ind(1) = 0
1355          vg(0) = vg_surface
1356          DO  k = 1, nzt+1
1357             IF ( i < 11 )  THEN
1358                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
1359                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
1360                   gradient = vg_vertical_gradient(i) / 100.0_wp
1361                   vg_vertical_gradient_level_ind(i) = k - 1
1362                   i = i + 1
1363                ENDIF
1364             ENDIF
1365             IF ( gradient /= 0.0_wp )  THEN
1366                IF ( k /= 1 )  THEN
1367                   vg(k) = vg(k-1) + dzu(k) * gradient
1368                ELSE
1369                   vg(k) = vg_surface + dzu(k) * gradient
1370                ENDIF
1371             ELSE
1372                vg(k) = vg(k-1)
1373             ENDIF
1374          ENDDO
1375
1376       ELSE
1377
1378          vg_vertical_gradient_level_ind(1) = nzt+1
1379          vg(nzt+1) = vg_surface
1380          DO  k = nzt, nzb, -1
1381             IF ( i < 11 )  THEN
1382                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.              &
1383                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
1384                   gradient = vg_vertical_gradient(i) / 100.0_wp
1385                   vg_vertical_gradient_level_ind(i) = k + 1
1386                   i = i + 1
1387                ENDIF
1388             ENDIF
1389             IF ( gradient /= 0.0_wp )  THEN
1390                IF ( k /= nzt )  THEN
1391                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1392                ELSE
1393                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1394                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1395                ENDIF
1396             ELSE
1397                vg(k) = vg(k+1)
1398             ENDIF
1399          ENDDO
1400
1401       ENDIF
1402
1403!
1404!--    In case of no given gradients for vg, choose a zero gradient
1405       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1406          vg_vertical_gradient_level(1) = 0.0_wp
1407       ENDIF
1408
1409!
1410!--    Let the initial wind profiles be the calculated ug/vg profiles or
1411!--    interpolate them from wind profile data (if given)
1412       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1413
1414          u_init = ug
1415          v_init = vg
1416
1417       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1418
1419          IF ( uv_heights(1) /= 0.0_wp )  THEN
1420             message_string = 'uv_heights(1) must be 0.0'
1421             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1422          ENDIF
1423
1424          IF ( omega /= 0.0_wp )  THEN
1425             message_string = 'Coriolis force must be switched off (by setting omega=0.0)' //  &
1426                              ' when prescribing the forcing by u_profile and v_profile'
1427             CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 )
1428          ENDIF
1429
1430          use_prescribed_profile_data = .TRUE.
1431
1432          kk = 1
1433          u_init(0) = 0.0_wp
1434          v_init(0) = 0.0_wp
1435
1436          DO  k = 1, nz+1
1437
1438             IF ( kk < 100 )  THEN
1439                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
1440                   kk = kk + 1
1441                   IF ( kk == 100 )  EXIT
1442                ENDDO
1443             ENDIF
1444
1445             IF ( kk < 100  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
1446                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1447                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1448                                       ( u_profile(kk+1) - u_profile(kk) )
1449                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1450                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1451                                       ( v_profile(kk+1) - v_profile(kk) )
1452             ELSE
1453                u_init(k) = u_profile(kk)
1454                v_init(k) = v_profile(kk)
1455             ENDIF
1456
1457          ENDDO
1458
1459       ELSE
1460
1461          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1462          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1463
1464       ENDIF
1465
1466!
1467!--    Compute initial temperature profile using the given temperature gradients
1468       IF (  .NOT.  neutral )  THEN
1469          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
1470                                       pt_vertical_gradient_level,              &
1471                                       pt_vertical_gradient, pt_init,           &
1472                                       pt_surface, bc_pt_t_val )
1473       ENDIF
1474!
1475!--    Compute initial humidity profile using the given humidity gradients
1476       IF ( humidity )  THEN
1477          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
1478                                       q_vertical_gradient_level,              &
1479                                       q_vertical_gradient, q_init,            &
1480                                       q_surface, bc_q_t_val )
1481       ENDIF
1482!
1483!--    Compute initial scalar profile using the given scalar gradients
1484       IF ( passive_scalar )  THEN
1485          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
1486                                       s_vertical_gradient_level,              &
1487                                       s_vertical_gradient, s_init,            &
1488                                       s_surface, bc_s_t_val )
1489       ENDIF
1490!
1491!--    If required, compute initial salinity profile using the given salinity
1492!--    gradients
1493       IF ( ocean )  THEN
1494          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
1495                                       sa_vertical_gradient_level,              &
1496                                       sa_vertical_gradient, sa_init,           &
1497                                       sa_surface, dum )
1498       ENDIF
1499
1500
1501    ENDIF
1502
1503!
1504!-- Check if the control parameter use_subsidence_tendencies is used correctly
1505    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
1506       message_string = 'The usage of use_subsidence_tendencies ' //           &
1507                            'requires large_scale_subsidence = .T..'
1508       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1509    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1510       message_string = 'The usage of use_subsidence_tendencies ' //           &
1511                            'requires large_scale_forcing = .T..'
1512       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1513    ENDIF
1514
1515!
1516!-- Initialize large scale subsidence if required
1517    If ( large_scale_subsidence )  THEN
1518       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
1519                                     .NOT.  large_scale_forcing )  THEN
1520          CALL init_w_subsidence
1521       ENDIF
1522!
1523!--    In case large_scale_forcing is used, profiles for subsidence velocity
1524!--    are read in from file LSF_DATA
1525
1526       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
1527            .NOT.  large_scale_forcing )  THEN
1528          message_string = 'There is no default large scale vertical ' //      &
1529                           'velocity profile set. Specify the subsidence ' //  &
1530                           'velocity profile via subs_vertical_gradient ' //   &
1531                           'and subs_vertical_gradient_level.'
1532          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1533       ENDIF
1534    ELSE
1535        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1536           message_string = 'Enable usage of large scale subsidence by ' //    &
1537                            'setting large_scale_subsidence = .T..'
1538          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1539        ENDIF
1540    ENDIF
1541
1542!
1543!-- Compute Coriolis parameter
1544    f  = 2.0_wp * omega * SIN( phi / 180.0_wp * pi )
1545    fs = 2.0_wp * omega * COS( phi / 180.0_wp * pi )
1546
1547!
1548!-- Check and set buoyancy related parameters and switches
1549    IF ( reference_state == 'horizontal_average' )  THEN
1550       CONTINUE
1551    ELSEIF ( reference_state == 'initial_profile' )  THEN
1552       use_initial_profile_as_reference = .TRUE.
1553    ELSEIF ( reference_state == 'single_value' )  THEN
1554       use_single_reference_value = .TRUE.
1555       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1556       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1557    ELSE
1558       message_string = 'illegal value for reference_state: "' //              &
1559                        TRIM( reference_state ) // '"'
1560       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1561    ENDIF
1562
1563!
1564!-- Sign of buoyancy/stability terms
1565    IF ( ocean )  atmos_ocean_sign = -1.0_wp
1566
1567!
1568!-- Ocean version must use flux boundary conditions at the top
1569    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1570       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
1571       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1572    ENDIF
1573
1574!
1575!-- In case of a given slope, compute the relevant quantities
1576    IF ( alpha_surface /= 0.0_wp )  THEN
1577       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1578          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
1579                                     ' ) must be < 90.0'
1580          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1581       ENDIF
1582       sloping_surface = .TRUE.
1583       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1584       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1585    ENDIF
1586
1587!
1588!-- Check time step and cfl_factor
1589    IF ( dt /= -1.0_wp )  THEN
1590       IF ( dt <= 0.0_wp  .AND.  dt /= -1.0_wp )  THEN
1591          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1592          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1593       ENDIF
1594       dt_3d = dt
1595       dt_fixed = .TRUE.
1596    ENDIF
1597
1598    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1599       IF ( cfl_factor == -1.0_wp )  THEN
1600          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1601             cfl_factor = 0.8_wp
1602          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1603             cfl_factor = 0.9_wp
1604          ELSE
1605             cfl_factor = 0.9_wp
1606          ENDIF
1607       ELSE
1608          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
1609                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
1610          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1611       ENDIF
1612    ENDIF
1613
1614!
1615!-- Store simulated time at begin
1616    simulated_time_at_begin = simulated_time
1617
1618!
1619!-- Store reference time for coupled runs and change the coupling flag,
1620!-- if ...
1621    IF ( simulated_time == 0.0_wp )  THEN
1622       IF ( coupling_start_time == 0.0_wp )  THEN
1623          time_since_reference_point = 0.0_wp
1624       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1625          run_coupled = .FALSE.
1626       ENDIF
1627    ENDIF
1628
1629!
1630!-- Set wind speed in the Galilei-transformed system
1631    IF ( galilei_transformation )  THEN
1632       IF ( use_ug_for_galilei_tr                    .AND.                     &
1633            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1634            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
1635            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1636            vg_vertical_gradient(1) == 0.0_wp )  THEN
1637          u_gtrans = ug_surface * 0.6_wp
1638          v_gtrans = vg_surface * 0.6_wp
1639       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1640                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1641                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1642          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
1643                           ' with galilei transformation'
1644          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1645       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1646                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1647                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1648          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
1649                           ' with galilei transformation'
1650          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1651       ELSE
1652          message_string = 'variable translation speed used for galilei-' //   &
1653             'transformation, which may cause & instabilities in stably ' //   &
1654             'stratified regions'
1655          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1656       ENDIF
1657    ENDIF
1658
1659!
1660!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1661!-- fluxes have to be used in the diffusion-terms
1662    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
1663
1664!
1665!-- Check boundary conditions and set internal variables:
1666!-- Attention: the lateral boundary conditions have been already checked in
1667!-- parin
1668!
1669!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1670!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1671!-- and tools do not work with non-cyclic boundary conditions.
1672    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1673       IF ( psolver(1:9) /= 'multigrid' )  THEN
1674          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1675                           'psolver = "' // TRIM( psolver ) // '"'
1676          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1677       ENDIF
1678       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
1679            momentum_advec /= 'ws-scheme' )  THEN
1680
1681          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1682                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1683          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1684       ENDIF
1685       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
1686            scalar_advec /= 'ws-scheme' )  THEN
1687          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1688                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1689          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1690       ENDIF
1691       IF ( galilei_transformation )  THEN
1692          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1693                           'galilei_transformation = .T.'
1694          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1695       ENDIF
1696    ENDIF
1697
1698!
1699!-- Bottom boundary condition for the turbulent Kinetic energy
1700    IF ( bc_e_b == 'neumann' )  THEN
1701       ibc_e_b = 1
1702    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1703       ibc_e_b = 2
1704       IF ( .NOT. constant_flux_layer )  THEN
1705          bc_e_b = 'neumann'
1706          ibc_e_b = 1
1707          message_string = 'boundary condition bc_e_b changed to "' //         &
1708                           TRIM( bc_e_b ) // '"'
1709          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1710       ENDIF
1711    ELSE
1712       message_string = 'unknown boundary condition: bc_e_b = "' //            &
1713                        TRIM( bc_e_b ) // '"'
1714       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1715    ENDIF
1716
1717!
1718!-- Boundary conditions for perturbation pressure
1719    IF ( bc_p_b == 'dirichlet' )  THEN
1720       ibc_p_b = 0
1721    ELSEIF ( bc_p_b == 'neumann' )  THEN
1722       ibc_p_b = 1
1723    ELSE
1724       message_string = 'unknown boundary condition: bc_p_b = "' //            &
1725                        TRIM( bc_p_b ) // '"'
1726       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1727    ENDIF
1728
1729    IF ( bc_p_t == 'dirichlet' )  THEN
1730       ibc_p_t = 0
1731!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
1732    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
1733       ibc_p_t = 1
1734    ELSE
1735       message_string = 'unknown boundary condition: bc_p_t = "' //            &
1736                        TRIM( bc_p_t ) // '"'
1737       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1738    ENDIF
1739
1740!
1741!-- Boundary conditions for potential temperature
1742    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1743       ibc_pt_b = 2
1744    ELSE
1745       IF ( bc_pt_b == 'dirichlet' )  THEN
1746          ibc_pt_b = 0
1747       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1748          ibc_pt_b = 1
1749       ELSE
1750          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
1751                           TRIM( bc_pt_b ) // '"'
1752          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1753       ENDIF
1754    ENDIF
1755
1756    IF ( bc_pt_t == 'dirichlet' )  THEN
1757       ibc_pt_t = 0
1758    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1759       ibc_pt_t = 1
1760    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1761       ibc_pt_t = 2
1762    ELSEIF ( bc_pt_t == 'nested' )  THEN
1763       ibc_pt_t = 3
1764    ELSE
1765       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
1766                        TRIM( bc_pt_t ) // '"'
1767       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1768    ENDIF
1769
1770    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.                        &
1771         surface_heatflux == 9999999.9_wp )  THEN
1772       message_string = 'wall_heatflux additionally requires ' //     &
1773                        'setting of surface_heatflux'
1774       CALL message( 'check_parameters', 'PA0443', 1, 2, 0, 6, 0 )
1775    ENDIF
1776
1777!
1778!   This IF clause needs revision, got too complex!!
1779    IF ( surface_heatflux == 9999999.9_wp  )  THEN
1780       constant_heatflux = .FALSE.
1781       IF ( large_scale_forcing  .OR.  land_surface  .OR.  urban_surface )  THEN
1782          IF ( ibc_pt_b == 0 )  THEN
1783             constant_heatflux = .FALSE.
1784          ELSEIF ( ibc_pt_b == 1 )  THEN
1785             constant_heatflux = .TRUE.
1786             IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.   &
1787                  .NOT.  land_surface )  THEN
1788                surface_heatflux = shf_surf(1)
1789             ELSE
1790                surface_heatflux = 0.0_wp
1791             ENDIF
1792          ENDIF
1793       ENDIF
1794    ELSE
1795        constant_heatflux = .TRUE.
1796        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.        &
1797             large_scale_forcing  .AND.  .NOT.  land_surface )  THEN
1798           surface_heatflux = shf_surf(1)
1799        ENDIF
1800    ENDIF
1801
1802    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
1803
1804    IF ( neutral )  THEN
1805
1806       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
1807            surface_heatflux /= 9999999.9_wp )  THEN
1808          message_string = 'heatflux must not be set for pure neutral flow'
1809          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1810       ENDIF
1811
1812       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
1813       THEN
1814          message_string = 'heatflux must not be set for pure neutral flow'
1815          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1816       ENDIF
1817
1818    ENDIF
1819
1820    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
1821         top_momentumflux_v /= 9999999.9_wp )  THEN
1822       constant_top_momentumflux = .TRUE.
1823    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
1824           top_momentumflux_v == 9999999.9_wp ) )  THEN
1825       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
1826                        'must be set'
1827       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1828    ENDIF
1829
1830!
1831!-- A given surface temperature implies Dirichlet boundary condition for
1832!-- temperature. In this case specification of a constant heat flux is
1833!-- forbidden.
1834    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
1835         surface_heatflux /= 0.0_wp )  THEN
1836       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1837                        '& is not allowed with constant_heatflux = .TRUE.'
1838       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1839    ENDIF
1840    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
1841       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
1842               'wed with pt_surface_initial_change (/=0) = ',                  &
1843               pt_surface_initial_change
1844       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1845    ENDIF
1846
1847!
1848!-- A given temperature at the top implies Dirichlet boundary condition for
1849!-- temperature. In this case specification of a constant heat flux is
1850!-- forbidden.
1851    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
1852         top_heatflux /= 0.0_wp )  THEN
1853       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1854                        '" is not allowed with constant_top_heatflux = .TRUE.'
1855       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1856    ENDIF
1857
1858!
1859!-- Boundary conditions for salinity
1860    IF ( ocean )  THEN
1861       IF ( bc_sa_t == 'dirichlet' )  THEN
1862          ibc_sa_t = 0
1863       ELSEIF ( bc_sa_t == 'neumann' )  THEN
1864          ibc_sa_t = 1
1865       ELSE
1866          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
1867                           TRIM( bc_sa_t ) // '"'
1868          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1869       ENDIF
1870
1871       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
1872       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
1873          message_string = 'boundary condition: bc_sa_t = "' //                &
1874                           TRIM( bc_sa_t ) // '" requires to set ' //          &
1875                           'top_salinityflux'
1876          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1877       ENDIF
1878
1879!
1880!--    A fixed salinity at the top implies Dirichlet boundary condition for
1881!--    salinity. In this case specification of a constant salinity flux is
1882!--    forbidden.
1883       IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.             &
1884            top_salinityflux /= 0.0_wp )  THEN
1885          message_string = 'boundary condition: bc_sa_t = "' //                &
1886                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
1887                           'top_salinityflux /= 0.0'
1888          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
1889       ENDIF
1890
1891    ENDIF
1892
1893!
1894!-- Set boundary conditions for total water content
1895    IF ( humidity )  THEN
1896
1897       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.                        &
1898            surface_waterflux == 9999999.9_wp )  THEN
1899          message_string = 'wall_humidityflux additionally requires ' //     &
1900                           'setting of surface_waterflux'
1901          CALL message( 'check_parameters', 'PA0444', 1, 2, 0, 6, 0 )
1902       ENDIF
1903
1904       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
1905                            'PA0071', 'PA0072' )
1906
1907       IF ( surface_waterflux == 9999999.9_wp  )  THEN
1908          constant_waterflux = .FALSE.
1909          IF ( large_scale_forcing .OR. land_surface )  THEN
1910             IF ( ibc_q_b == 0 )  THEN
1911                constant_waterflux = .FALSE.
1912             ELSEIF ( ibc_q_b == 1 )  THEN
1913                constant_waterflux = .TRUE.
1914                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1915                   surface_waterflux = qsws_surf(1)
1916                ENDIF
1917             ENDIF
1918          ENDIF
1919       ELSE
1920          constant_waterflux = .TRUE.
1921          IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.      &
1922               large_scale_forcing )  THEN
1923             surface_waterflux = qsws_surf(1)
1924          ENDIF
1925       ENDIF
1926
1927       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074',        &
1928                              constant_waterflux, q_surface_initial_change )
1929
1930    ENDIF
1931
1932    IF ( passive_scalar )  THEN
1933
1934       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.                        &
1935            surface_scalarflux == 9999999.9_wp )  THEN
1936          message_string = 'wall_scalarflux additionally requires ' //     &
1937                           'setting of surface_scalarflux'
1938          CALL message( 'check_parameters', 'PA0445', 1, 2, 0, 6, 0 )
1939       ENDIF
1940
1941       IF ( surface_scalarflux == 9999999.9_wp )  constant_scalarflux = .FALSE.
1942
1943       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,             &
1944                            'PA0071', 'PA0072' )
1945
1946       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074',        &
1947                              constant_scalarflux, s_surface_initial_change )
1948
1949       IF ( top_scalarflux == 9999999.9_wp )  constant_top_scalarflux = .FALSE.
1950!
1951!--    A fixed scalar concentration at the top implies Dirichlet boundary
1952!--    condition for scalar. Hence, in this case specification of a constant
1953!--    scalar flux is forbidden.
1954       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
1955               .AND.  top_scalarflux /= 0.0_wp )  THEN
1956          message_string = 'boundary condition: bc_s_t = "' //                 &
1957                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
1958                           'top_scalarflux /= 0.0'
1959          CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
1960       ENDIF
1961    ENDIF
1962!
1963!-- Boundary conditions for horizontal components of wind speed
1964    IF ( bc_uv_b == 'dirichlet' )  THEN
1965       ibc_uv_b = 0
1966    ELSEIF ( bc_uv_b == 'neumann' )  THEN
1967       ibc_uv_b = 1
1968       IF ( constant_flux_layer )  THEN
1969          message_string = 'boundary condition: bc_uv_b = "' //                &
1970               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
1971               // ' = .TRUE.'
1972          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
1973       ENDIF
1974    ELSE
1975       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
1976                        TRIM( bc_uv_b ) // '"'
1977       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
1978    ENDIF
1979!
1980!-- In case of coupled simulations u and v at the ground in atmosphere will be
1981!-- assigned with the u and v values of the ocean surface
1982    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1983       ibc_uv_b = 2
1984    ENDIF
1985
1986    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1987       bc_uv_t = 'neumann'
1988       ibc_uv_t = 1
1989    ELSE
1990       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
1991          ibc_uv_t = 0
1992          IF ( bc_uv_t == 'dirichlet_0' )  THEN
1993!
1994!--          Velocities for the initial u,v-profiles are set zero at the top
1995!--          in case of dirichlet_0 conditions
1996             u_init(nzt+1)    = 0.0_wp
1997             v_init(nzt+1)    = 0.0_wp
1998          ENDIF
1999       ELSEIF ( bc_uv_t == 'neumann' )  THEN
2000          ibc_uv_t = 1
2001       ELSEIF ( bc_uv_t == 'nested' )  THEN
2002          ibc_uv_t = 3
2003       ELSE
2004          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
2005                           TRIM( bc_uv_t ) // '"'
2006          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
2007       ENDIF
2008    ENDIF
2009
2010!
2011!-- Compute and check, respectively, the Rayleigh Damping parameter
2012    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
2013       rayleigh_damping_factor = 0.0_wp
2014    ELSE
2015       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                             &
2016            rayleigh_damping_factor > 1.0_wp )  THEN
2017          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
2018                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
2019          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
2020       ENDIF
2021    ENDIF
2022
2023    IF ( rayleigh_damping_height == -1.0_wp )  THEN
2024       IF (  .NOT.  ocean )  THEN
2025          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
2026       ELSE
2027          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
2028       ENDIF
2029    ELSE
2030       IF (  .NOT.  ocean )  THEN
2031          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
2032               rayleigh_damping_height > zu(nzt) )  THEN
2033             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2034                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
2035             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2036          ENDIF
2037       ELSE
2038          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
2039               rayleigh_damping_height < zu(nzb) )  THEN
2040             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2041                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
2042             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2043          ENDIF
2044       ENDIF
2045    ENDIF
2046
2047!
2048!-- Check number of chosen statistic regions
2049    IF ( statistic_regions < 0 )  THEN
2050       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
2051                   statistic_regions+1, ' is not allowed'
2052       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
2053    ENDIF
2054    IF ( normalizing_region > statistic_regions  .OR.                          &
2055         normalizing_region < 0)  THEN
2056       WRITE ( message_string, * ) 'normalizing_region = ',                    &
2057                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
2058                ' (value of statistic_regions)'
2059       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
2060    ENDIF
2061
2062!
2063!-- Set the default intervals for data output, if necessary
2064!-- NOTE: dt_dosp has already been set in spectra_parin
2065    IF ( dt_data_output /= 9999999.9_wp )  THEN
2066       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
2067       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
2068       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
2069       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
2070       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
2071       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
2072       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
2073       DO  mid = 1, max_masks
2074          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
2075       ENDDO
2076    ENDIF
2077
2078!
2079!-- Set the default skip time intervals for data output, if necessary
2080    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
2081                                       skip_time_dopr    = skip_time_data_output
2082    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
2083                                       skip_time_do2d_xy = skip_time_data_output
2084    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
2085                                       skip_time_do2d_xz = skip_time_data_output
2086    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
2087                                       skip_time_do2d_yz = skip_time_data_output
2088    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
2089                                       skip_time_do3d    = skip_time_data_output
2090    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
2091                                skip_time_data_output_av = skip_time_data_output
2092    DO  mid = 1, max_masks
2093       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
2094                                skip_time_domask(mid)    = skip_time_data_output
2095    ENDDO
2096
2097!
2098!-- Check the average intervals (first for 3d-data, then for profiles)
2099    IF ( averaging_interval > dt_data_output_av )  THEN
2100       WRITE( message_string, * )  'averaging_interval = ',                    &
2101             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
2102       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
2103    ENDIF
2104
2105    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
2106       averaging_interval_pr = averaging_interval
2107    ENDIF
2108
2109    IF ( averaging_interval_pr > dt_dopr )  THEN
2110       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
2111             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
2112       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
2113    ENDIF
2114
2115!
2116!-- Set the default interval for profiles entering the temporal average
2117    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
2118       dt_averaging_input_pr = dt_averaging_input
2119    ENDIF
2120
2121!
2122!-- Set the default interval for the output of timeseries to a reasonable
2123!-- value (tries to minimize the number of calls of flow_statistics)
2124    IF ( dt_dots == 9999999.9_wp )  THEN
2125       IF ( averaging_interval_pr == 0.0_wp )  THEN
2126          dt_dots = MIN( dt_run_control, dt_dopr )
2127       ELSE
2128          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
2129       ENDIF
2130    ENDIF
2131
2132!
2133!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
2134    IF ( dt_averaging_input > averaging_interval )  THEN
2135       WRITE( message_string, * )  'dt_averaging_input = ',                    &
2136                dt_averaging_input, ' must be <= averaging_interval = ',       &
2137                averaging_interval
2138       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
2139    ENDIF
2140
2141    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
2142       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
2143                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
2144                averaging_interval_pr
2145       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2146    ENDIF
2147
2148!
2149!-- Set the default value for the integration interval of precipitation amount
2150    IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
2151       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
2152          precipitation_amount_interval = dt_do2d_xy
2153       ELSE
2154          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
2155             WRITE( message_string, * )  'precipitation_amount_interval = ',   &
2156                 precipitation_amount_interval, ' must not be larger than ',   &
2157                 'dt_do2d_xy = ', dt_do2d_xy
2158             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
2159          ENDIF
2160       ENDIF
2161    ENDIF
2162
2163!
2164!-- Determine the number of output profiles and check whether they are
2165!-- permissible
2166    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2167
2168       dopr_n = dopr_n + 1
2169       i = dopr_n
2170
2171!
2172!--    Determine internal profile number (for hom, homs)
2173!--    and store height levels
2174       SELECT CASE ( TRIM( data_output_pr(i) ) )
2175
2176          CASE ( 'u', '#u' )
2177             dopr_index(i) = 1
2178             dopr_unit(i)  = 'm/s'
2179             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2180             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2181                dopr_initial_index(i) = 5
2182                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2183                data_output_pr(i)     = data_output_pr(i)(2:)
2184             ENDIF
2185
2186          CASE ( 'v', '#v' )
2187             dopr_index(i) = 2
2188             dopr_unit(i)  = 'm/s'
2189             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2190             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2191                dopr_initial_index(i) = 6
2192                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2193                data_output_pr(i)     = data_output_pr(i)(2:)
2194             ENDIF
2195
2196          CASE ( 'w' )
2197             dopr_index(i) = 3
2198             dopr_unit(i)  = 'm/s'
2199             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2200
2201          CASE ( 'pt', '#pt' )
2202             IF ( .NOT. cloud_physics ) THEN
2203                dopr_index(i) = 4
2204                dopr_unit(i)  = 'K'
2205                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2206                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2207                   dopr_initial_index(i) = 7
2208                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2209                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2210                   data_output_pr(i)     = data_output_pr(i)(2:)
2211                ENDIF
2212             ELSE
2213                dopr_index(i) = 43
2214                dopr_unit(i)  = 'K'
2215                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2216                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2217                   dopr_initial_index(i) = 28
2218                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2219                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2220                   data_output_pr(i)     = data_output_pr(i)(2:)
2221                ENDIF
2222             ENDIF
2223
2224          CASE ( 'e' )
2225             dopr_index(i)  = 8
2226             dopr_unit(i)   = 'm2/s2'
2227             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2228             hom(nzb,2,8,:) = 0.0_wp
2229
2230          CASE ( 'km', '#km' )
2231             dopr_index(i)  = 9
2232             dopr_unit(i)   = 'm2/s'
2233             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2234             hom(nzb,2,9,:) = 0.0_wp
2235             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2236                dopr_initial_index(i) = 23
2237                hom(:,2,23,:)         = hom(:,2,9,:)
2238                data_output_pr(i)     = data_output_pr(i)(2:)
2239             ENDIF
2240
2241          CASE ( 'kh', '#kh' )
2242             dopr_index(i)   = 10
2243             dopr_unit(i)    = 'm2/s'
2244             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2245             hom(nzb,2,10,:) = 0.0_wp
2246             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2247                dopr_initial_index(i) = 24
2248                hom(:,2,24,:)         = hom(:,2,10,:)
2249                data_output_pr(i)     = data_output_pr(i)(2:)
2250             ENDIF
2251
2252          CASE ( 'l', '#l' )
2253             dopr_index(i)   = 11
2254             dopr_unit(i)    = 'm'
2255             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2256             hom(nzb,2,11,:) = 0.0_wp
2257             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2258                dopr_initial_index(i) = 25
2259                hom(:,2,25,:)         = hom(:,2,11,:)
2260                data_output_pr(i)     = data_output_pr(i)(2:)
2261             ENDIF
2262
2263          CASE ( 'w"u"' )
2264             dopr_index(i) = 12
2265             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2266             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2267             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
2268
2269          CASE ( 'w*u*' )
2270             dopr_index(i) = 13
2271             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2272             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2273
2274          CASE ( 'w"v"' )
2275             dopr_index(i) = 14
2276             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2277             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2278             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
2279
2280          CASE ( 'w*v*' )
2281             dopr_index(i) = 15
2282             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2283             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2284
2285          CASE ( 'w"pt"' )
2286             dopr_index(i) = 16
2287             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2288             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2289
2290          CASE ( 'w*pt*' )
2291             dopr_index(i) = 17
2292             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2293             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2294
2295          CASE ( 'wpt' )
2296             dopr_index(i) = 18
2297             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2298             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2299
2300          CASE ( 'wu' )
2301             dopr_index(i) = 19
2302             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2303             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2304             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
2305
2306          CASE ( 'wv' )
2307             dopr_index(i) = 20
2308             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2309             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2310             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
2311
2312          CASE ( 'w*pt*BC' )
2313             dopr_index(i) = 21
2314             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2315             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2316
2317          CASE ( 'wptBC' )
2318             dopr_index(i) = 22
2319             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2320             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2321
2322          CASE ( 'sa', '#sa' )
2323             IF ( .NOT. ocean )  THEN
2324                message_string = 'data_output_pr = ' // &
2325                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2326                                 'lemented for ocean = .FALSE.'
2327                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2328             ELSE
2329                dopr_index(i) = 23
2330                dopr_unit(i)  = 'psu'
2331                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2332                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2333                   dopr_initial_index(i) = 26
2334                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2335                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2336                   data_output_pr(i)     = data_output_pr(i)(2:)
2337                ENDIF
2338             ENDIF
2339
2340          CASE ( 'u*2' )
2341             dopr_index(i) = 30
2342             dopr_unit(i)  = 'm2/s2'
2343             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2344
2345          CASE ( 'v*2' )
2346             dopr_index(i) = 31
2347             dopr_unit(i)  = 'm2/s2'
2348             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2349
2350          CASE ( 'w*2' )
2351             dopr_index(i) = 32
2352             dopr_unit(i)  = 'm2/s2'
2353             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2354
2355          CASE ( 'pt*2' )
2356             dopr_index(i) = 33
2357             dopr_unit(i)  = 'K2'
2358             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2359
2360          CASE ( 'e*' )
2361             dopr_index(i) = 34
2362             dopr_unit(i)  = 'm2/s2'
2363             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2364
2365          CASE ( 'w*2pt*' )
2366             dopr_index(i) = 35
2367             dopr_unit(i)  = 'K m2/s2'
2368             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2369
2370          CASE ( 'w*pt*2' )
2371             dopr_index(i) = 36
2372             dopr_unit(i)  = 'K2 m/s'
2373             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2374
2375          CASE ( 'w*e*' )
2376             dopr_index(i) = 37
2377             dopr_unit(i)  = 'm3/s3'
2378             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2379
2380          CASE ( 'w*3' )
2381             dopr_index(i) = 38
2382             dopr_unit(i)  = 'm3/s3'
2383             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2384
2385          CASE ( 'Sw' )
2386             dopr_index(i) = 39
2387             dopr_unit(i)  = 'none'
2388             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2389
2390          CASE ( 'p' )
2391             dopr_index(i) = 40
2392             dopr_unit(i)  = 'Pa'
2393             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2394
2395          CASE ( 'q', '#q' )
2396             IF ( .NOT. humidity )  THEN
2397                message_string = 'data_output_pr = ' // &
2398                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2399                                 'lemented for humidity = .FALSE.'
2400                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2401             ELSE
2402                dopr_index(i) = 41
2403                dopr_unit(i)  = 'kg/kg'
2404                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2405                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2406                   dopr_initial_index(i) = 26
2407                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2408                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2409                   data_output_pr(i)     = data_output_pr(i)(2:)
2410                ENDIF
2411             ENDIF
2412
2413          CASE ( 's', '#s' )
2414             IF ( .NOT. passive_scalar )  THEN
2415                message_string = 'data_output_pr = ' //                        &
2416                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2417                                 'lemented for passive_scalar = .FALSE.'
2418                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2419             ELSE
2420                dopr_index(i) = 117
2421                dopr_unit(i)  = 'kg/m3'
2422                hom(:,2,117,:) = SPREAD( zu, 2, statistic_regions+1 )
2423                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2424                   dopr_initial_index(i) = 117
2425                   hom(:,2,117,:)        = SPREAD( zu, 2, statistic_regions+1 )
2426                   hom(nzb,2,117,:)      = 0.0_wp    ! because zu(nzb) is negative
2427                   data_output_pr(i)     = data_output_pr(i)(2:)
2428                ENDIF
2429             ENDIF
2430
2431          CASE ( 'qv', '#qv' )
2432             IF ( .NOT. cloud_physics ) THEN
2433                dopr_index(i) = 41
2434                dopr_unit(i)  = 'kg/kg'
2435                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2436                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2437                   dopr_initial_index(i) = 26
2438                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2439                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2440                   data_output_pr(i)     = data_output_pr(i)(2:)
2441                ENDIF
2442             ELSE
2443                dopr_index(i) = 42
2444                dopr_unit(i)  = 'kg/kg'
2445                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2446                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2447                   dopr_initial_index(i) = 27
2448                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2449                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2450                   data_output_pr(i)     = data_output_pr(i)(2:)
2451                ENDIF
2452             ENDIF
2453
2454          CASE ( 'lpt', '#lpt' )
2455             IF ( .NOT. cloud_physics ) THEN
2456                message_string = 'data_output_pr = ' //                        &
2457                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2458                                 'lemented for cloud_physics = .FALSE.'
2459                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2460             ELSE
2461                dopr_index(i) = 4
2462                dopr_unit(i)  = 'K'
2463                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2464                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2465                   dopr_initial_index(i) = 7
2466                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2467                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2468                   data_output_pr(i)     = data_output_pr(i)(2:)
2469                ENDIF
2470             ENDIF
2471
2472          CASE ( 'vpt', '#vpt' )
2473             dopr_index(i) = 44
2474             dopr_unit(i)  = 'K'
2475             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2476             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2477                dopr_initial_index(i) = 29
2478                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2479                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2480                data_output_pr(i)     = data_output_pr(i)(2:)
2481             ENDIF
2482
2483          CASE ( 'w"vpt"' )
2484             dopr_index(i) = 45
2485             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2486             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2487
2488          CASE ( 'w*vpt*' )
2489             dopr_index(i) = 46
2490             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2491             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2492
2493          CASE ( 'wvpt' )
2494             dopr_index(i) = 47
2495             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2496             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2497
2498          CASE ( 'w"q"' )
2499             IF (  .NOT.  humidity )  THEN
2500                message_string = 'data_output_pr = ' //                        &
2501                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2502                                 'lemented for humidity = .FALSE.'
2503                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2504             ELSE
2505                dopr_index(i) = 48
2506                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2507                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2508             ENDIF
2509
2510          CASE ( 'w*q*' )
2511             IF (  .NOT.  humidity )  THEN
2512                message_string = 'data_output_pr = ' //                        &
2513                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2514                                 'lemented for humidity = .FALSE.'
2515                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2516             ELSE
2517                dopr_index(i) = 49
2518                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2519                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2520             ENDIF
2521
2522          CASE ( 'wq' )
2523             IF (  .NOT.  humidity )  THEN
2524                message_string = 'data_output_pr = ' //                        &
2525                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2526                                 'lemented for humidity = .FALSE.'
2527                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2528             ELSE
2529                dopr_index(i) = 50
2530                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2531                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2532             ENDIF
2533
2534          CASE ( 'w"s"' )
2535             IF (  .NOT.  passive_scalar )  THEN
2536                message_string = 'data_output_pr = ' //                        &
2537                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2538                                 'lemented for passive_scalar = .FALSE.'
2539                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2540             ELSE
2541                dopr_index(i) = 119
2542                dopr_unit(i)  = 'kg/m3 m/s'
2543                hom(:,2,119,:) = SPREAD( zw, 2, statistic_regions+1 )
2544             ENDIF
2545
2546          CASE ( 'w*s*' )
2547             IF (  .NOT.  passive_scalar )  THEN
2548                message_string = 'data_output_pr = ' //                        &
2549                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2550                                 'lemented for passive_scalar = .FALSE.'
2551                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2552             ELSE
2553                dopr_index(i) = 116
2554                dopr_unit(i)  = 'kg/m3 m/s'
2555                hom(:,2,116,:) = SPREAD( zw, 2, statistic_regions+1 )
2556             ENDIF
2557
2558          CASE ( 'ws' )
2559             IF (  .NOT.  passive_scalar )  THEN
2560                message_string = 'data_output_pr = ' //                        &
2561                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2562                                 'lemented for passive_scalar = .FALSE.'
2563                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2564             ELSE
2565                dopr_index(i) = 120
2566                dopr_unit(i)  = 'kg/m3 m/s'
2567                hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
2568             ENDIF
2569
2570          CASE ( 'w"qv"' )
2571             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2572                dopr_index(i) = 48
2573                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2574                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2575             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2576                dopr_index(i) = 51
2577                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2578                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2579             ELSE
2580                message_string = 'data_output_pr = ' //                        &
2581                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2582                                 'lemented for cloud_physics = .FALSE. an&' // &
2583                                 'd humidity = .FALSE.'
2584                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2585             ENDIF
2586
2587          CASE ( 'w*qv*' )
2588             IF ( humidity  .AND.  .NOT. cloud_physics )                       &
2589             THEN
2590                dopr_index(i) = 49
2591                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2592                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2593             ELSEIF( humidity .AND. cloud_physics ) THEN
2594                dopr_index(i) = 52
2595                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2596                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2597             ELSE
2598                message_string = 'data_output_pr = ' //                        &
2599                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2600                                 'lemented for cloud_physics = .FALSE. an&' // &
2601                                 'd humidity = .FALSE.'
2602                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2603             ENDIF
2604
2605          CASE ( 'wqv' )
2606             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2607                dopr_index(i) = 50
2608                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2609                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2610             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2611                dopr_index(i) = 53
2612                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2613                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2614             ELSE
2615                message_string = 'data_output_pr = ' //                        &
2616                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2617                                 'lemented for cloud_physics = .FALSE. an&' // &
2618                                 'd humidity = .FALSE.'
2619                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2620             ENDIF
2621
2622          CASE ( 'ql' )
2623             IF (  .NOT.  cloud_physics  .AND.  .NOT.  cloud_droplets )  THEN
2624                message_string = 'data_output_pr = ' //                        &
2625                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2626                                 'lemented for cloud_physics = .FALSE. or'  // &
2627                                 '&cloud_droplets = .FALSE.'
2628                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2629             ELSE
2630                dopr_index(i) = 54
2631                dopr_unit(i)  = 'kg/kg'
2632                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2633             ENDIF
2634
2635          CASE ( 'w*u*u*:dz' )
2636             dopr_index(i) = 55
2637             dopr_unit(i)  = 'm2/s3'
2638             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2639
2640          CASE ( 'w*p*:dz' )
2641             dopr_index(i) = 56
2642             dopr_unit(i)  = 'm2/s3'
2643             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2644
2645          CASE ( 'w"e:dz' )
2646             dopr_index(i) = 57
2647             dopr_unit(i)  = 'm2/s3'
2648             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2649
2650
2651          CASE ( 'u"pt"' )
2652             dopr_index(i) = 58
2653             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2654             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2655
2656          CASE ( 'u*pt*' )
2657             dopr_index(i) = 59
2658             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2659             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2660
2661          CASE ( 'upt_t' )
2662             dopr_index(i) = 60
2663             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2664             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2665
2666          CASE ( 'v"pt"' )
2667             dopr_index(i) = 61
2668             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2669             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2670
2671          CASE ( 'v*pt*' )
2672             dopr_index(i) = 62
2673             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2674             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2675
2676          CASE ( 'vpt_t' )
2677             dopr_index(i) = 63
2678             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2679             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2680
2681          CASE ( 'rho_ocean' )
2682             IF (  .NOT.  ocean ) THEN
2683                message_string = 'data_output_pr = ' //                        &
2684                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2685                                 'lemented for ocean = .FALSE.'
2686                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2687             ELSE
2688                dopr_index(i) = 64
2689                dopr_unit(i)  = 'kg/m3'
2690                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2691                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2692                   dopr_initial_index(i) = 77
2693                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2694                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
2695                   data_output_pr(i)     = data_output_pr(i)(2:)
2696                ENDIF
2697             ENDIF
2698
2699          CASE ( 'w"sa"' )
2700             IF (  .NOT.  ocean ) THEN
2701                message_string = 'data_output_pr = ' //                        &
2702                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2703                                 'lemented for ocean = .FALSE.'
2704                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2705             ELSE
2706                dopr_index(i) = 65
2707                dopr_unit(i)  = 'psu m/s'
2708                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2709             ENDIF
2710
2711          CASE ( 'w*sa*' )
2712             IF (  .NOT. ocean  ) THEN
2713                message_string = 'data_output_pr = ' //                        &
2714                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2715                                 'lemented for ocean = .FALSE.'
2716                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2717             ELSE
2718                dopr_index(i) = 66
2719                dopr_unit(i)  = 'psu m/s'
2720                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2721             ENDIF
2722
2723          CASE ( 'wsa' )
2724             IF (  .NOT.  ocean ) THEN
2725                message_string = 'data_output_pr = ' //                        &
2726                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2727                                 'lemented for ocean = .FALSE.'
2728                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2729             ELSE
2730                dopr_index(i) = 67
2731                dopr_unit(i)  = 'psu m/s'
2732                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2733             ENDIF
2734
2735          CASE ( 'w*p*' )
2736             dopr_index(i) = 68
2737             dopr_unit(i)  = 'm3/s3'
2738             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2739
2740          CASE ( 'w"e' )
2741             dopr_index(i) = 69
2742             dopr_unit(i)  = 'm3/s3'
2743             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2744
2745          CASE ( 'q*2' )
2746             IF (  .NOT.  humidity )  THEN
2747                message_string = 'data_output_pr = ' //                        &
2748                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2749                                 'lemented for humidity = .FALSE.'
2750                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2751             ELSE
2752                dopr_index(i) = 70
2753                dopr_unit(i)  = 'kg2/kg2'
2754                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2755             ENDIF
2756
2757          CASE ( 'prho' )
2758             IF (  .NOT.  ocean ) THEN
2759                message_string = 'data_output_pr = ' //                        &
2760                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2761                                 'lemented for ocean = .FALSE.'
2762                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2763             ELSE
2764                dopr_index(i) = 71
2765                dopr_unit(i)  = 'kg/m3'
2766                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2767             ENDIF
2768
2769          CASE ( 'hyp' )
2770             dopr_index(i) = 72
2771             dopr_unit(i)  = 'hPa'
2772             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2773
2774          CASE ( 'rho_air' )
2775             dopr_index(i)  = 119
2776             dopr_unit(i)   = 'kg/m3'
2777             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
2778
2779          CASE ( 'rho_air_zw' )
2780             dopr_index(i)  = 120
2781             dopr_unit(i)   = 'kg/m3'
2782             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
2783
2784          CASE ( 'nc' )
2785             IF (  .NOT.  cloud_physics )  THEN
2786                message_string = 'data_output_pr = ' //                        &
2787                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2788                                 'lemented for cloud_physics = .FALSE.'
2789                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2790             ELSEIF ( .NOT.  microphysics_morrison )  THEN
2791                message_string = 'data_output_pr = ' //                        &
2792                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2793                                 'lemented for cloud_scheme /= morrison'
2794                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2795             ELSE
2796                dopr_index(i) = 89
2797                dopr_unit(i)  = '1/m3'
2798                hom(:,2,89,:)  = SPREAD( zu, 2, statistic_regions+1 )
2799             ENDIF
2800
2801          CASE ( 'nr' )
2802             IF (  .NOT.  cloud_physics )  THEN
2803                message_string = 'data_output_pr = ' //                        &
2804                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2805                                 'lemented for cloud_physics = .FALSE.'
2806                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2807             ELSEIF ( .NOT.  microphysics_seifert )  THEN
2808                message_string = 'data_output_pr = ' //                        &
2809                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2810                                 'lemented for cloud_scheme /= seifert_beheng'
2811                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2812             ELSE
2813                dopr_index(i) = 73
2814                dopr_unit(i)  = '1/m3'
2815                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
2816             ENDIF
2817
2818          CASE ( 'qr' )
2819             IF (  .NOT.  cloud_physics )  THEN
2820                message_string = 'data_output_pr = ' //                        &
2821                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2822                                 'lemented for cloud_physics = .FALSE.'
2823                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2824             ELSEIF ( .NOT.  microphysics_seifert )  THEN
2825                message_string = 'data_output_pr = ' //                        &
2826                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2827                                 'lemented for cloud_scheme /= seifert_beheng'
2828                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2829             ELSE
2830                dopr_index(i) = 74
2831                dopr_unit(i)  = 'kg/kg'
2832                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
2833             ENDIF
2834
2835          CASE ( 'qc' )
2836             IF (  .NOT.  cloud_physics )  THEN
2837                message_string = 'data_output_pr = ' //                        &
2838                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2839                                 'lemented for cloud_physics = .FALSE.'
2840                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2841             ELSE
2842                dopr_index(i) = 75
2843                dopr_unit(i)  = 'kg/kg'
2844                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
2845             ENDIF
2846
2847          CASE ( 'prr' )
2848             IF (  .NOT.  cloud_physics )  THEN
2849                message_string = 'data_output_pr = ' //                        &
2850                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2851                                 'lemented for cloud_physics = .FALSE.'
2852                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2853             ELSEIF ( microphysics_sat_adjust )  THEN
2854                message_string = 'data_output_pr = ' //                        &
2855                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
2856                                 'ilable for cloud_scheme = saturation_adjust'
2857                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
2858             ELSE
2859                dopr_index(i) = 76
2860                dopr_unit(i)  = 'kg/kg m/s'
2861                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
2862             ENDIF
2863
2864          CASE ( 'ug' )
2865             dopr_index(i) = 78
2866             dopr_unit(i)  = 'm/s'
2867             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2868
2869          CASE ( 'vg' )
2870             dopr_index(i) = 79
2871             dopr_unit(i)  = 'm/s'
2872             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2873
2874          CASE ( 'w_subs' )
2875             IF (  .NOT.  large_scale_subsidence )  THEN
2876                message_string = 'data_output_pr = ' //                        &
2877                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2878                                 'lemented for large_scale_subsidence = .FALSE.'
2879                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
2880             ELSE
2881                dopr_index(i) = 80
2882                dopr_unit(i)  = 'm/s'
2883                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
2884             ENDIF
2885
2886          CASE ( 's*2' )
2887             IF (  .NOT.  passive_scalar )  THEN
2888                message_string = 'data_output_pr = ' //                        &
2889                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2890                                 'lemented for passive_scalar = .FALSE.'
2891                CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
2892             ELSE
2893                dopr_index(i) = 118
2894                dopr_unit(i)  = 'kg2/kg2'
2895                hom(:,2,118,:) = SPREAD( zu, 2, statistic_regions+1 )
2896             ENDIF
2897
2898
2899
2900          CASE DEFAULT
2901
2902             CALL lsm_check_data_output_pr( data_output_pr(i), i, unit,        &
2903                                            dopr_unit(i) )
2904
2905             CALL lsf_nudging_check_data_output_pr( data_output_pr(i), i,      &
2906                                                    unit, dopr_unit(i) )
2907
2908             IF ( unit == 'illegal' )  THEN
2909                CALL radiation_check_data_output_pr( data_output_pr(i), i,     &
2910                                                     unit, dopr_unit(i) )
2911             ENDIF
2912
2913             IF ( unit == 'illegal' )  THEN
2914                unit = ''
2915                CALL user_check_data_output_pr( data_output_pr(i), i, unit )
2916             ENDIF
2917
2918             IF ( unit == 'illegal' )  THEN
2919                IF ( data_output_pr_user(1) /= ' ' )  THEN
2920                   message_string = 'illegal value for data_output_pr or ' //  &
2921                                    'data_output_pr_user = "' //               &
2922                                    TRIM( data_output_pr(i) ) // '"'
2923                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2924                ELSE
2925                   message_string = 'illegal value for data_output_pr = "' //  &
2926                                    TRIM( data_output_pr(i) ) // '"'
2927                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2928                ENDIF
2929             ENDIF
2930
2931       END SELECT
2932
2933    ENDDO
2934
2935
2936!
2937!-- Append user-defined data output variables to the standard data output
2938    IF ( data_output_user(1) /= ' ' )  THEN
2939       i = 1
2940       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2941          i = i + 1
2942       ENDDO
2943       j = 1
2944       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
2945          IF ( i > 500 )  THEN
2946             message_string = 'number of output quantitities given by data' // &
2947                '_output and data_output_user exceeds the limit of 500'
2948             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2949          ENDIF
2950          data_output(i) = data_output_user(j)
2951          i = i + 1
2952          j = j + 1
2953       ENDDO
2954    ENDIF
2955
2956!
2957!-- Check and set steering parameters for 2d/3d data output and averaging
2958    i   = 1
2959    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2960!
2961!--    Check for data averaging
2962       ilen = LEN_TRIM( data_output(i) )
2963       j = 0                                                 ! no data averaging
2964       IF ( ilen > 3 )  THEN
2965          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
2966             j = 1                                           ! data averaging
2967             data_output(i) = data_output(i)(1:ilen-3)
2968          ENDIF
2969       ENDIF
2970!
2971!--    Check for cross section or volume data
2972       ilen = LEN_TRIM( data_output(i) )
2973       k = 0                                                   ! 3d data
2974       var = data_output(i)(1:ilen)
2975       IF ( ilen > 3 )  THEN
2976          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
2977               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
2978               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2979             k = 1                                             ! 2d data
2980             var = data_output(i)(1:ilen-3)
2981          ENDIF
2982       ENDIF
2983
2984!
2985!--    Check for allowed value and set units
2986       SELECT CASE ( TRIM( var ) )
2987
2988          CASE ( 'e' )
2989             IF ( constant_diffusion )  THEN
2990                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2991                                 'res constant_diffusion = .FALSE.'
2992                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
2993             ENDIF
2994             unit = 'm2/s2'
2995
2996          CASE ( 'lpt' )
2997             IF (  .NOT.  cloud_physics )  THEN
2998                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
2999                         'res cloud_physics = .TRUE.'
3000                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3001             ENDIF
3002             unit = 'K'
3003
3004          CASE ( 'nc' )
3005             IF (  .NOT.  cloud_physics )  THEN
3006                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3007                         'res cloud_physics = .TRUE.'
3008                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3009             ELSEIF ( .NOT.  microphysics_morrison )  THEN
3010                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3011                         'res = microphysics morrison '
3012                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3013             ENDIF
3014             unit = '1/m3'
3015
3016          CASE ( 'nr' )
3017             IF (  .NOT.  cloud_physics )  THEN
3018                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3019                         'res cloud_physics = .TRUE.'
3020                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3021             ELSEIF ( .NOT.  microphysics_seifert )  THEN
3022                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3023                         'res cloud_scheme = seifert_beheng'
3024                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3025             ENDIF
3026             unit = '1/m3'
3027
3028          CASE ( 'pc', 'pr' )
3029             IF (  .NOT.  particle_advection )  THEN
3030                message_string = 'output of "' // TRIM( var ) // '" requir' // &
3031                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
3032                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
3033             ENDIF
3034             IF ( TRIM( var ) == 'pc' )  unit = 'number'
3035             IF ( TRIM( var ) == 'pr' )  unit = 'm'
3036
3037          CASE ( 'prr' )
3038             IF (  .NOT.  cloud_physics )  THEN
3039                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3040                         'res cloud_physics = .TRUE.'
3041                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3042             ELSEIF ( microphysics_sat_adjust )  THEN
3043                message_string = 'output of "' // TRIM( var ) // '" is ' //  &
3044                         'not available for cloud_scheme = saturation_adjust'
3045                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
3046             ENDIF
3047             unit = 'kg/kg m/s'
3048
3049          CASE ( 'q', 'vpt' )
3050             IF (  .NOT.  humidity )  THEN
3051                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3052                                 'res humidity = .TRUE.'
3053                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
3054             ENDIF
3055             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
3056             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
3057
3058          CASE ( 'qc' )
3059             IF (  .NOT.  cloud_physics )  THEN
3060                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3061                         'res cloud_physics = .TRUE.'
3062                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3063             ENDIF
3064             unit = 'kg/kg'
3065
3066          CASE ( 'ql' )
3067             IF ( .NOT.  ( cloud_physics  .OR.  cloud_droplets ) )  THEN
3068                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3069                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
3070                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
3071             ENDIF
3072             unit = 'kg/kg'
3073
3074          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
3075             IF (  .NOT.  cloud_droplets )  THEN
3076                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3077                                 'res cloud_droplets = .TRUE.'
3078                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
3079             ENDIF
3080             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
3081             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
3082             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
3083
3084          CASE ( 'qr' )
3085             IF (  .NOT.  cloud_physics )  THEN
3086                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3087                         'res cloud_physics = .TRUE.'
3088                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3089             ELSEIF ( .NOT.  microphysics_seifert ) THEN
3090                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3091                         'res cloud_scheme = seifert_beheng'
3092                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3093             ENDIF
3094             unit = 'kg/kg'
3095
3096          CASE ( 'qv' )
3097             IF (  .NOT.  cloud_physics )  THEN
3098                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3099                                 'res cloud_physics = .TRUE.'
3100                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3101             ENDIF
3102             unit = 'kg/kg'
3103
3104          CASE ( 'rho_ocean' )
3105             IF (  .NOT.  ocean )  THEN
3106                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3107                                 'res ocean = .TRUE.'
3108                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3109             ENDIF
3110             unit = 'kg/m3'
3111
3112          CASE ( 's' )
3113             IF (  .NOT.  passive_scalar )  THEN
3114                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3115                                 'res passive_scalar = .TRUE.'
3116                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3117             ENDIF
3118             unit = 'conc'
3119
3120          CASE ( 'sa' )
3121             IF (  .NOT.  ocean )  THEN
3122                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3123                                 'res ocean = .TRUE.'
3124                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3125             ENDIF
3126             unit = 'psu'
3127
3128          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3129             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3130             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3131             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3132             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3133             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3134             CONTINUE
3135
3136          CASE ( 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'shf*', 'ssws*', 't*', &
3137                 'u*', 'z0*', 'z0h*', 'z0q*' )
3138             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3139                message_string = 'illegal value for data_output: "' //         &
3140                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3141                                 'cross sections are allowed for this value'
3142                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3143             ENDIF
3144
3145             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
3146                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3147                                 'res cloud_physics = .TRUE.'
3148                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3149             ENDIF
3150             IF ( TRIM( var ) == 'pra*'  .AND.                                 &
3151                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
3152                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3153                                 'res cloud_scheme = kessler or seifert_beheng'
3154                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3155             ENDIF
3156             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
3157                message_string = 'temporal averaging of precipitation ' //     &
3158                          'amount "' // TRIM( var ) // '" is not possible'
3159                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
3160             ENDIF
3161             IF ( TRIM( var ) == 'prr*'  .AND.                                 &
3162                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
3163                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3164                                 'res cloud_scheme = kessler or seifert_beheng'
3165                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3166             ENDIF
3167             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
3168                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3169                                 'res humidity = .TRUE.'
3170                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3171             ENDIF
3172             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
3173                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3174                                 'res passive_scalar = .TRUE.'
3175                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
3176             ENDIF
3177
3178             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
3179             IF ( TRIM( var ) == 'ol*'   )   unit = 'm'
3180             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
3181             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
3182             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3183             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3184             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
3185             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3186             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3187             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3188             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3189
3190          CASE DEFAULT
3191
3192             CALL lsm_check_data_output ( var, unit, i, ilen, k )
3193
3194             IF ( unit == 'illegal' )  THEN
3195                CALL radiation_check_data_output( var, unit, i, ilen, k )
3196             ENDIF
3197
3198!
3199!--          Block of urban surface model outputs
3200             IF ( unit == 'illegal' .AND. urban_surface .AND. var(1:4) == 'usm_' ) THEN
3201                 CALL usm_check_data_output( var, unit )
3202             ENDIF
3203
3204!
3205!--          Block of plant canopy model outputs
3206             IF ( unit == 'illegal' .AND. plant_canopy .AND. var(1:4) == 'pcm_' ) THEN
3207                 CALL pcm_check_data_output( var, unit )
3208             ENDIF
3209
3210             IF ( unit == 'illegal' )  THEN
3211                unit = ''
3212                CALL user_check_data_output( var, unit )
3213             ENDIF
3214
3215             IF ( unit == 'illegal' )  THEN
3216                IF ( data_output_user(1) /= ' ' )  THEN
3217                   message_string = 'illegal value for data_output or ' //     &
3218                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3219                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3220                ELSE
3221                   message_string = 'illegal value for data_output = "' //     &
3222                                    TRIM( data_output(i) ) // '"'
3223                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3224                ENDIF
3225             ENDIF
3226
3227       END SELECT
3228!
3229!--    Set the internal steering parameters appropriately
3230       IF ( k == 0 )  THEN
3231          do3d_no(j)              = do3d_no(j) + 1
3232          do3d(j,do3d_no(j))      = data_output(i)
3233          do3d_unit(j,do3d_no(j)) = unit
3234       ELSE
3235          do2d_no(j)              = do2d_no(j) + 1
3236          do2d(j,do2d_no(j))      = data_output(i)
3237          do2d_unit(j,do2d_no(j)) = unit
3238          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3239             data_output_xy(j) = .TRUE.
3240          ENDIF
3241          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3242             data_output_xz(j) = .TRUE.
3243          ENDIF
3244          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3245             data_output_yz(j) = .TRUE.
3246          ENDIF
3247       ENDIF
3248
3249       IF ( j == 1 )  THEN
3250!
3251!--       Check, if variable is already subject to averaging
3252          found = .FALSE.
3253          DO  k = 1, doav_n
3254             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3255          ENDDO
3256
3257          IF ( .NOT. found )  THEN
3258             doav_n = doav_n + 1
3259             doav(doav_n) = var
3260          ENDIF
3261       ENDIF
3262
3263       i = i + 1
3264    ENDDO
3265
3266!
3267!-- Averaged 2d or 3d output requires that an averaging interval has been set
3268    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3269       WRITE( message_string, * )  'output of averaged quantity "',            &
3270                                   TRIM( doav(1) ), '_av" requires to set a ', &
3271                                   'non-zero & averaging interval'
3272       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3273    ENDIF
3274
3275!
3276!-- Check sectional planes and store them in one shared array
3277    IF ( ANY( section_xy > nz + 1 ) )  THEN
3278       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3279       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3280    ENDIF
3281    IF ( ANY( section_xz > ny + 1 ) )  THEN
3282       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3283       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3284    ENDIF
3285    IF ( ANY( section_yz > nx + 1 ) )  THEN
3286       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3287       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3288    ENDIF
3289    section(:,1) = section_xy
3290    section(:,2) = section_xz
3291    section(:,3) = section_yz
3292
3293!
3294!-- Upper plot limit for 2D vertical sections
3295    IF ( z_max_do2d == -1.0_wp )  z_max_do2d = zu(nzt)
3296    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
3297       WRITE( message_string, * )  'z_max_do2d = ', z_max_do2d,                &
3298                    ' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
3299                    ' (zu(nzt))'
3300       CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
3301    ENDIF
3302
3303!
3304!-- Upper plot limit for 3D arrays
3305    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3306
3307!
3308!-- Set output format string (used in header)
3309    SELECT CASE ( netcdf_data_format )
3310       CASE ( 1 )
3311          netcdf_data_format_string = 'netCDF classic'
3312       CASE ( 2 )
3313          netcdf_data_format_string = 'netCDF 64bit offset'
3314       CASE ( 3 )
3315          netcdf_data_format_string = 'netCDF4/HDF5'
3316       CASE ( 4 )
3317          netcdf_data_format_string = 'netCDF4/HDF5 classic'
3318       CASE ( 5 )
3319          netcdf_data_format_string = 'parallel netCDF4/HDF5'
3320       CASE ( 6 )
3321          netcdf_data_format_string = 'parallel netCDF4/HDF5 classic'
3322
3323    END SELECT
3324
3325!
3326!-- Check mask conditions
3327    DO mid = 1, max_masks
3328       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
3329            data_output_masks_user(mid,1) /= ' ' ) THEN
3330          masks = masks + 1
3331       ENDIF
3332    ENDDO
3333
3334    IF ( masks < 0  .OR.  masks > max_masks )  THEN
3335       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3336            '<= ', max_masks, ' (=max_masks)'
3337       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3338    ENDIF
3339    IF ( masks > 0 )  THEN
3340       mask_scale(1) = mask_scale_x
3341       mask_scale(2) = mask_scale_y
3342       mask_scale(3) = mask_scale_z
3343       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3344          WRITE( message_string, * )                                           &
3345               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3346               'must be > 0.0'
3347          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3348       ENDIF
3349!
3350!--    Generate masks for masked data output
3351!--    Parallel netcdf output is not tested so far for masked data, hence
3352!--    netcdf_data_format is switched back to non-paralell output.
3353       netcdf_data_format_save = netcdf_data_format
3354       IF ( netcdf_data_format > 4 )  THEN
3355          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3356          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3357          message_string = 'netCDF file formats '//                            &
3358                           '5 (parallel netCDF 4) and ' //                     &
3359                           '6 (parallel netCDF 4 Classic model) '//            &
3360                           '&are currently not supported (not yet tested) ' // &
3361                           'for masked data.&Using respective non-parallel' // &
3362                           ' output for masked data.'
3363          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3364       ENDIF
3365       CALL init_masks
3366       netcdf_data_format = netcdf_data_format_save
3367    ENDIF
3368
3369!
3370!-- Check the NetCDF data format
3371    IF ( netcdf_data_format > 2 )  THEN
3372#if defined( __netcdf4 )
3373       CONTINUE
3374#else
3375       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3376                        'cpp-directive __netcdf4 given & switch '  //          &
3377                        'back to 64-bit offset format'
3378       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3379       netcdf_data_format = 2
3380#endif
3381    ENDIF
3382    IF ( netcdf_data_format > 4 )  THEN
3383#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3384       CONTINUE
3385#else
3386       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3387                        'cpp-directive __netcdf4_parallel given & switch '  // &
3388                        'back to netCDF4 non-parallel output'
3389       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3390       netcdf_data_format = netcdf_data_format - 2
3391#endif
3392    ENDIF
3393
3394!
3395!-- Calculate fixed number of output time levels for parallel netcdf output.
3396!-- The time dimension has to be defined as limited for parallel output,
3397!-- because otherwise the I/O performance drops significantly.
3398    IF ( netcdf_data_format > 4 )  THEN
3399
3400!
3401!--    Check if any of the follwoing data output interval is 0.0s, which is
3402!--    not allowed for parallel output.
3403       CALL check_dt_do( dt_do3d,    'dt_do3d'    )
3404       CALL check_dt_do( dt_do2d_xy, 'dt_do2d_xy' )
3405       CALL check_dt_do( dt_do2d_xz, 'dt_do2d_xz' )
3406       CALL check_dt_do( dt_do2d_yz, 'dt_do2d_yz' )
3407
3408       ntdim_3d(0)    = INT( ( end_time - skip_time_do3d ) / dt_do3d )
3409       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3410       ntdim_3d(1)    = INT( ( end_time - skip_time_data_output_av )           &
3411                             / dt_data_output_av )
3412       ntdim_2d_xy(0) = INT( ( end_time - skip_time_do2d_xy ) / dt_do2d_xy )
3413       ntdim_2d_xz(0) = INT( ( end_time - skip_time_do2d_xz ) / dt_do2d_xz )
3414       ntdim_2d_yz(0) = INT( ( end_time - skip_time_do2d_yz ) / dt_do2d_yz )
3415       IF ( do2d_at_begin )  THEN
3416          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3417          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3418          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3419       ENDIF
3420       ntdim_2d_xy(1) = ntdim_3d(1)
3421       ntdim_2d_xz(1) = ntdim_3d(1)
3422       ntdim_2d_yz(1) = ntdim_3d(1)
3423
3424    ENDIF
3425
3426!
3427!-- Check, whether a constant diffusion coefficient shall be used
3428    IF ( km_constant /= -1.0_wp )  THEN
3429       IF ( km_constant < 0.0_wp )  THEN
3430          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3431          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3432       ELSE
3433          IF ( prandtl_number < 0.0_wp )  THEN
3434             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3435                                         ' < 0.0'
3436             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3437          ENDIF
3438          constant_diffusion = .TRUE.
3439
3440          IF ( constant_flux_layer )  THEN
3441             message_string = 'constant_flux_layer is not allowed with fixed ' &
3442                              // 'value of km'
3443             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3444          ENDIF
3445       ENDIF
3446    ENDIF
3447
3448!
3449!-- In case of non-cyclic lateral boundaries and a damping layer for the
3450!-- potential temperature, check the width of the damping layer
3451    IF ( bc_lr /= 'cyclic' ) THEN
3452       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3453            pt_damping_width > REAL( (nx+1) * dx ) )  THEN
3454          message_string = 'pt_damping_width out of range'
3455          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3456       ENDIF
3457    ENDIF
3458
3459    IF ( bc_ns /= 'cyclic' )  THEN
3460       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3461            pt_damping_width > REAL( (ny+1) * dy ) )  THEN
3462          message_string = 'pt_damping_width out of range'
3463          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3464       ENDIF
3465    ENDIF
3466
3467!
3468!-- Check value range for zeta = z/L
3469    IF ( zeta_min >= zeta_max )  THEN
3470       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
3471                                   'than zeta_max = ', zeta_max
3472       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3473    ENDIF
3474
3475!
3476!-- Check random generator
3477    IF ( (random_generator /= 'system-specific'      .AND.                     &
3478          random_generator /= 'random-parallel'   )  .AND.                     &
3479          random_generator /= 'numerical-recipes' )  THEN
3480       message_string = 'unknown random generator: random_generator = "' //    &
3481                        TRIM( random_generator ) // '"'
3482       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3483    ENDIF
3484
3485!
3486!-- Determine upper and lower hight level indices for random perturbations
3487    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3488       IF ( ocean )  THEN
3489          disturbance_level_b     = zu((nzt*2)/3)
3490          disturbance_level_ind_b = ( nzt * 2 ) / 3
3491       ELSE
3492          disturbance_level_b     = zu(nzb+3)
3493          disturbance_level_ind_b = nzb + 3
3494       ENDIF
3495    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3496       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3497                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3498       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3499    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3500       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3501                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3502       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3503    ELSE
3504       DO  k = 3, nzt-2
3505          IF ( disturbance_level_b <= zu(k) )  THEN
3506             disturbance_level_ind_b = k
3507             EXIT
3508          ENDIF
3509       ENDDO
3510    ENDIF
3511
3512    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3513       IF ( ocean )  THEN
3514          disturbance_level_t     = zu(nzt-3)
3515          disturbance_level_ind_t = nzt - 3
3516       ELSE
3517          disturbance_level_t     = zu(nzt/3)
3518          disturbance_level_ind_t = nzt / 3
3519       ENDIF
3520    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3521       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3522                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3523       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3524    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3525       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3526                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3527                   disturbance_level_b
3528       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3529    ELSE
3530       DO  k = 3, nzt-2
3531          IF ( disturbance_level_t <= zu(k) )  THEN
3532             disturbance_level_ind_t = k
3533             EXIT
3534          ENDIF
3535       ENDDO
3536    ENDIF
3537
3538!
3539!-- Check again whether the levels determined this way are ok.
3540!-- Error may occur at automatic determination and too few grid points in
3541!-- z-direction.
3542    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3543       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3544                disturbance_level_ind_t, ' must be >= disturbance_level_ind_b = ', &
3545                disturbance_level_ind_b
3546       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3547    ENDIF
3548
3549!
3550!-- Determine the horizontal index range for random perturbations.
3551!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3552!-- near the inflow and the perturbation area is further limited to ...(1)
3553!-- after the initial phase of the flow.
3554
3555    IF ( bc_lr /= 'cyclic' )  THEN
3556       IF ( inflow_disturbance_begin == -1 )  THEN
3557          inflow_disturbance_begin = MIN( 10, nx/2 )
3558       ENDIF
3559       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3560       THEN
3561          message_string = 'inflow_disturbance_begin out of range'
3562          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3563       ENDIF
3564       IF ( inflow_disturbance_end == -1 )  THEN
3565          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3566       ENDIF
3567       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3568       THEN
3569          message_string = 'inflow_disturbance_end out of range'
3570          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3571       ENDIF
3572    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3573       IF ( inflow_disturbance_begin == -1 )  THEN
3574          inflow_disturbance_begin = MIN( 10, ny/2 )
3575       ENDIF
3576       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3577       THEN
3578          message_string = 'inflow_disturbance_begin out of range'
3579          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3580       ENDIF
3581       IF ( inflow_disturbance_end == -1 )  THEN
3582          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3583       ENDIF
3584       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3585       THEN
3586          message_string = 'inflow_disturbance_end out of range'
3587          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3588       ENDIF
3589    ENDIF
3590
3591    IF ( random_generator == 'random-parallel' )  THEN
3592       dist_nxl = nxl;  dist_nxr = nxr
3593       dist_nys = nys;  dist_nyn = nyn
3594       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3595          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3596          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
3597       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3598          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3599          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
3600       ELSEIF ( bc_lr == 'nested' )  THEN
3601          dist_nxl = MAX( inflow_disturbance_begin, nxl )
3602       ENDIF
3603       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3604          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3605          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
3606       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3607          dist_nys    = MAX( inflow_disturbance_begin, nys )
3608          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
3609       ELSEIF ( bc_ns == 'nested' )  THEN
3610          dist_nys = MAX( inflow_disturbance_begin, nys )
3611       ENDIF
3612    ELSE
3613       dist_nxl = 0;  dist_nxr = nx
3614       dist_nys = 0;  dist_nyn = ny
3615       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3616          dist_nxr    = nx - inflow_disturbance_begin
3617          dist_nxl(1) = nx - inflow_disturbance_end
3618       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3619          dist_nxl    = inflow_disturbance_begin
3620          dist_nxr(1) = inflow_disturbance_end
3621       ENDIF
3622       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3623          dist_nyn    = ny - inflow_disturbance_begin
3624          dist_nys(1) = ny - inflow_disturbance_end
3625       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3626          dist_nys    = inflow_disturbance_begin
3627          dist_nyn(1) = inflow_disturbance_end
3628       ENDIF
3629    ENDIF
3630
3631!
3632!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3633!-- boundary (so far, a turbulent inflow is realized from the left side only)
3634    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3635       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
3636                        'condition at the inflow boundary'
3637       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3638    ENDIF
3639
3640!
3641!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3642!-- data from prerun in the first main run
3643    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3644         initializing_actions /= 'read_restart_data' )  THEN
3645       message_string = 'turbulent_inflow = .T. requires ' //                  &
3646                        'initializing_actions = ''cyclic_fill'' or ' //        &
3647                        'initializing_actions = ''read_restart_data'' '
3648       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3649    ENDIF
3650
3651!
3652!-- In case of turbulent inflow calculate the index of the recycling plane
3653    IF ( turbulent_inflow )  THEN
3654       IF ( recycling_width <= dx  .OR.  recycling_width >= nx * dx )  THEN
3655          WRITE( message_string, * )  'illegal value for recycling_width:', &
3656                                      ' ', recycling_width
3657          CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
3658       ENDIF
3659!
3660!--    Calculate the index
3661       recycling_plane = recycling_width / dx
3662!
3663!--    Because the y-shift is done with a distance of INT( npey / 2 ) no shift
3664!--    is possible if there is only one PE in y direction.
3665       IF ( recycling_yshift .AND. pdims(2) < 2 )  THEN
3666          WRITE( message_string, * )  'recycling_yshift = .T. requires more',  &
3667                                      ' than one processor in y direction'
3668          CALL message( 'check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
3669       ENDIF
3670    ENDIF
3671
3672
3673    IF ( turbulent_outflow )  THEN
3674!
3675!--    Turbulent outflow requires Dirichlet conditions at the respective inflow
3676!--    boundary (so far, a turbulent outflow is realized at the right side only)
3677       IF ( bc_lr /= 'dirichlet/radiation' )  THEN
3678          message_string = 'turbulent_outflow = .T. requires ' //              &
3679                           'bc_lr = "dirichlet/radiation"'
3680          CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
3681       ENDIF
3682!
3683!--    The ouflow-source plane must lay inside the model domain
3684       IF ( outflow_source_plane < dx  .OR.  &
3685            outflow_source_plane > nx * dx )  THEN
3686          WRITE( message_string, * )  'illegal value for outflow_source'//     &
3687                                      '_plane: ', outflow_source_plane
3688          CALL message( 'check_parameters', 'PA0145', 1, 2, 0, 6, 0 )
3689       ENDIF
3690    ENDIF
3691
3692!
3693!-- Determine damping level index for 1D model
3694    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
3695       IF ( damp_level_1d == -1.0_wp )  THEN
3696          damp_level_1d     = zu(nzt+1)
3697          damp_level_ind_1d = nzt + 1
3698       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
3699          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
3700                 ' must be >= 0.0 and <= ', zu(nzt+1), '(zu(nzt+1))'
3701          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
3702       ELSE
3703          DO  k = 1, nzt+1
3704             IF ( damp_level_1d <= zu(k) )  THEN
3705                damp_level_ind_1d = k
3706                EXIT
3707             ENDIF
3708          ENDDO
3709       ENDIF
3710    ENDIF
3711
3712!
3713!-- Check some other 1d-model parameters
3714    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
3715         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
3716       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
3717                        '" is unknown'
3718       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
3719    ENDIF
3720    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
3721         TRIM( dissipation_1d ) /= 'detering' )  THEN
3722       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
3723                        '" is unknown'
3724       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
3725    ENDIF
3726
3727!
3728!-- Set time for the next user defined restart (time_restart is the
3729!-- internal parameter for steering restart events)
3730    IF ( restart_time /= 9999999.9_wp )  THEN
3731       IF ( restart_time > time_since_reference_point )  THEN
3732          time_restart = restart_time
3733       ENDIF
3734    ELSE
3735!
3736!--    In case of a restart run, set internal parameter to default (no restart)
3737!--    if the NAMELIST-parameter restart_time is omitted
3738       time_restart = 9999999.9_wp
3739    ENDIF
3740
3741!
3742!-- Check pressure gradient conditions
3743    IF ( dp_external  .AND.  conserve_volume_flow )  THEN
3744       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
3745            'w are .TRUE. but one of them must be .FALSE.'
3746       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
3747    ENDIF
3748    IF ( dp_external )  THEN
3749       IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
3750          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
3751               ' of range'
3752          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
3753       ENDIF
3754       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
3755          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
3756               'ro, i.e. the external pressure gradient & will not be applied'
3757          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
3758       ENDIF
3759    ENDIF
3760    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
3761       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
3762            '.FALSE., i.e. the external pressure gradient & will not be applied'
3763       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
3764    ENDIF
3765    IF ( conserve_volume_flow )  THEN
3766       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
3767
3768          conserve_volume_flow_mode = 'initial_profiles'
3769
3770       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
3771            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
3772          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
3773               conserve_volume_flow_mode
3774          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
3775       ENDIF
3776       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
3777          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
3778          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
3779               'require  conserve_volume_flow_mode = ''initial_profiles'''
3780          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
3781       ENDIF
3782    ENDIF
3783    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
3784         ( .NOT. conserve_volume_flow  .OR.                                    &
3785         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
3786       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
3787            'conserve_volume_flow = .T. and ',                                 &
3788            'conserve_volume_flow_mode = ''bulk_velocity'''
3789       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
3790    ENDIF
3791
3792!
3793!-- Check particle attributes
3794    IF ( particle_color /= 'none' )  THEN
3795       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
3796            particle_color /= 'z' )  THEN
3797          message_string = 'illegal value for parameter particle_color: ' //   &
3798                           TRIM( particle_color)
3799          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
3800       ELSE
3801          IF ( color_interval(2) <= color_interval(1) )  THEN
3802             message_string = 'color_interval(2) <= color_interval(1)'
3803             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
3804          ENDIF
3805       ENDIF
3806    ENDIF
3807
3808    IF ( particle_dvrpsize /= 'none' )  THEN
3809       IF ( particle_dvrpsize /= 'absw' )  THEN
3810          message_string = 'illegal value for parameter particle_dvrpsize:' // &
3811                           ' ' // TRIM( particle_color)
3812          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
3813       ELSE
3814          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
3815             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
3816             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
3817          ENDIF
3818       ENDIF
3819    ENDIF
3820
3821!
3822!-- Prevent empty time records in volume, cross-section and masked data in case
3823!-- of non-parallel netcdf-output in restart runs
3824    IF ( netcdf_data_format < 5 )  THEN
3825       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
3826          do3d_time_count    = 0
3827          do2d_xy_time_count = 0
3828          do2d_xz_time_count = 0
3829          do2d_yz_time_count = 0
3830          domask_time_count  = 0
3831       ENDIF
3832    ENDIF
3833
3834!
3835!-- Check for valid setting of most_method
3836    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
3837         TRIM( most_method ) /= 'newton'    .AND.                              &
3838         TRIM( most_method ) /= 'lookup' )  THEN
3839       message_string = 'most_method = "' // TRIM( most_method ) //            &
3840                        '" is unknown'
3841       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
3842    ENDIF
3843
3844!
3845!-- Check roughness length, which has to be smaller than dz/2
3846    IF ( ( constant_flux_layer .OR.  &
3847           INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) &
3848         .AND. roughness_length >= 0.5 * dz )  THEN
3849       message_string = 'roughness_length must be smaller than dz/2'
3850       CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
3851    ENDIF
3852
3853    CALL location_message( 'finished', .TRUE. )
3854!
3855!-- Check &userpar parameters
3856    CALL user_check_parameters
3857
3858 CONTAINS
3859
3860!------------------------------------------------------------------------------!
3861! Description:
3862! ------------
3863!> Check the length of data output intervals. In case of parallel NetCDF output
3864!> the time levels of the output files need to be fixed. Therefore setting the
3865!> output interval to 0.0s (usually used to output each timestep) is not
3866!> possible as long as a non-fixed timestep is used.
3867!------------------------------------------------------------------------------!
3868
3869    SUBROUTINE check_dt_do( dt_do, dt_do_name )
3870
3871       IMPLICIT NONE
3872
3873       CHARACTER (LEN=*), INTENT (IN) :: dt_do_name !< parin variable name
3874
3875       REAL(wp), INTENT (INOUT)       :: dt_do      !< data output interval
3876
3877       IF ( dt_do == 0.0_wp )  THEN
3878          IF ( dt_fixed )  THEN
3879             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
3880                    'timestep is wanted (' // dt_do_name // ' = 0.0).&'//      &
3881                    'Setting the output interval to the fixed timestep '//     &
3882                    'dt = ', dt, 's.'
3883             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
3884             dt_do = dt
3885          ELSE
3886             message_string = dt_do_name // ' = 0.0 while using a ' //         &
3887                              'variable timestep and parallel netCDF4 ' //     &
3888                              'is not allowed.'
3889             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
3890          ENDIF
3891       ENDIF
3892
3893    END SUBROUTINE check_dt_do
3894
3895
3896
3897
3898!------------------------------------------------------------------------------!
3899! Description:
3900! ------------
3901!> Inititalizes the vertical profiles of scalar quantities.
3902!------------------------------------------------------------------------------!
3903
3904    SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,            &
3905                                       vertical_gradient_level,                &
3906                                       vertical_gradient,                      &
3907                                       pr_init, surface_value, bc_t_val )
3908
3909
3910       IMPLICIT NONE
3911
3912       INTEGER(iwp) ::  i     !< counter
3913       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
3914
3915       REAL(wp)     ::  bc_t_val      !< model top gradient
3916       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
3917       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
3918
3919
3920       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
3921       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
3922       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
3923
3924       i = 1
3925       gradient = 0.0_wp
3926
3927       IF (  .NOT.  ocean )  THEN
3928
3929          vertical_gradient_level_ind(1) = 0
3930          DO  k = 1, nzt+1
3931             IF ( i < 11 )  THEN
3932                IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
3933                     vertical_gradient_level(i) >= 0.0_wp )  THEN
3934                   gradient = vertical_gradient(i) / 100.0_wp
3935                   vertical_gradient_level_ind(i) = k - 1
3936                   i = i + 1
3937                ENDIF
3938             ENDIF
3939             IF ( gradient /= 0.0_wp )  THEN
3940                IF ( k /= 1 )  THEN
3941                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
3942                ELSE
3943                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
3944                ENDIF
3945             ELSE
3946                pr_init(k) = pr_init(k-1)
3947             ENDIF
3948   !
3949   !--       Avoid negative values
3950             IF ( pr_init(k) < 0.0_wp )  THEN
3951                pr_init(k) = 0.0_wp
3952             ENDIF
3953          ENDDO
3954
3955       ELSE
3956
3957          vertical_gradient_level_ind(1) = nzt+1
3958          DO  k = nzt, 0, -1
3959             IF ( i < 11 )  THEN
3960                IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
3961                     vertical_gradient_level(i) <= 0.0_wp )  THEN
3962                   gradient = vertical_gradient(i) / 100.0_wp
3963                   vertical_gradient_level_ind(i) = k + 1
3964                   i = i + 1
3965                ENDIF
3966             ENDIF
3967             IF ( gradient /= 0.0_wp )  THEN
3968                IF ( k /= nzt )  THEN
3969                   pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient
3970                ELSE
3971                   pr_init(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
3972                   pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
3973                ENDIF
3974             ELSE
3975                pr_init(k) = pr_init(k+1)
3976             ENDIF
3977   !
3978   !--       Avoid negative humidities
3979             IF ( pr_init(k) < 0.0_wp )  THEN
3980                pr_init(k) = 0.0_wp
3981             ENDIF
3982          ENDDO
3983
3984       ENDIF
3985
3986!
3987!--    In case of no given gradients, choose zero gradient conditions
3988       IF ( vertical_gradient_level(1) == -1.0_wp )  THEN
3989          vertical_gradient_level(1) = 0.0_wp
3990       ENDIF
3991!
3992!--    Store gradient at the top boundary for possible Neumann boundary condition
3993       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
3994
3995    END SUBROUTINE init_vertical_profiles
3996
3997
3998
3999!------------------------------------------------------------------------------!
4000! Description:
4001! ------------
4002!> Set the bottom and top boundary conditions for humidity and scalars.
4003!------------------------------------------------------------------------------!
4004
4005    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
4006
4007
4008       IMPLICIT NONE
4009
4010       CHARACTER (LEN=1)   ::  sq                       !<
4011       CHARACTER (LEN=*)   ::  bc_b
4012       CHARACTER (LEN=*)   ::  bc_t
4013       CHARACTER (LEN=*)   ::  err_nr_b
4014       CHARACTER (LEN=*)   ::  err_nr_t
4015
4016       INTEGER(iwp)        ::  ibc_b
4017       INTEGER(iwp)        ::  ibc_t
4018
4019!
4020!--    Set Integer flags and check for possilbe errorneous settings for bottom
4021!--    boundary condition
4022       IF ( bc_b == 'dirichlet' )  THEN
4023          ibc_b = 0
4024       ELSEIF ( bc_b == 'neumann' )  THEN
4025          ibc_b = 1
4026       ELSE
4027          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4028                           '_b ="' // TRIM( bc_b ) // '"'
4029          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
4030       ENDIF
4031!
4032!--    Set Integer flags and check for possilbe errorneous settings for top
4033!--    boundary condition
4034       IF ( bc_t == 'dirichlet' )  THEN
4035          ibc_t = 0
4036       ELSEIF ( bc_t == 'neumann' )  THEN
4037          ibc_t = 1
4038       ELSEIF ( bc_t == 'initial_gradient' )  THEN
4039          ibc_t = 2
4040       ELSEIF ( bc_t == 'nested' )  THEN
4041          ibc_t = 3
4042       ELSE
4043          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4044                           '_t ="' // TRIM( bc_t ) // '"'
4045          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
4046       ENDIF
4047
4048
4049    END SUBROUTINE set_bc_scalars
4050
4051
4052
4053!------------------------------------------------------------------------------!
4054! Description:
4055! ------------
4056!> Check for consistent settings of bottom boundary conditions for humidity
4057!> and scalars.
4058!------------------------------------------------------------------------------!
4059
4060    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b,                 &
4061                                 err_nr_1, err_nr_2,       &
4062                                 constant_flux, surface_initial_change )
4063
4064
4065       IMPLICIT NONE
4066
4067       CHARACTER (LEN=1)   ::  sq                       !<
4068       CHARACTER (LEN=*)   ::  bc_b
4069       CHARACTER (LEN=*)   ::  err_nr_1
4070       CHARACTER (LEN=*)   ::  err_nr_2
4071
4072       INTEGER(iwp)        ::  ibc_b
4073
4074       LOGICAL             ::  constant_flux
4075
4076       REAL(wp)            ::  surface_initial_change
4077
4078!
4079!--    A given surface value implies Dirichlet boundary condition for
4080!--    the respective quantity. In this case specification of a constant flux is
4081!--    forbidden. However, an exception is made for large-scale forcing as well
4082!--    as land-surface model.
4083       IF ( .NOT. land_surface  .AND.  .NOT. large_scale_forcing )  THEN
4084          IF ( ibc_b == 0  .AND.  constant_flux )  THEN
4085             message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' //  &
4086                              '= "' // TRIM( bc_b ) // '" is not allowed with ' // &
4087                              'prescribed surface flux'
4088             CALL message( 'check_parameters', err_nr_1, 1, 2, 0, 6, 0 )
4089          ENDIF
4090       ENDIF
4091       IF ( constant_flux  .AND.  surface_initial_change /= 0.0_wp )  THEN
4092          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
4093                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
4094                 surface_initial_change
4095          CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
4096       ENDIF
4097
4098
4099    END SUBROUTINE check_bc_scalars
4100
4101
4102
4103 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.