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

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

Merge of branch palm4u into trunk

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