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

Last change on this file since 2564 was 2564, checked in by Giersch, 7 years ago

Bugfix: Variable wind_turbine changed to module control_parameters

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