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

Last change on this file since 2088 was 2088, checked in by suehring, 5 years ago

Bugfixes in initial salinity profile and generic topography definition in case of ocean simulations

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