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

Last change on this file since 2249 was 2249, checked in by sward, 8 years ago

last commit documented

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