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

Last change on this file since 2766 was 2766, checked in by kanani, 6 years ago

Removal of chem directive, plus minor changes

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