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

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

Bugfix in the calculation of fixed number of output time levels for parallel netcdf output, error message related to reading sky view factors revised

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