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

Last change on this file since 2251 was 2251, checked in by Giersch, 4 years ago

last commit documented

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