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

Last change on this file since 2298 was 2292, checked in by schwenkel, 7 years ago

implementation of new bulk microphysics scheme

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