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

Last change on this file since 2765 was 2765, checked in by maronga, 4 years ago

major bugfix in calculation of aerodynamic resistance for vertical surface elements

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