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

Last change on this file since 2052 was 2052, checked in by gronemeier, 7 years ago

Bugfix: remove setting of default value for recycling_width

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