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

Last change on this file since 2200 was 2200, checked in by suehring, 7 years ago

Bugfixes in radiation_model; monotonic limiter in 5th order advection scheme removed

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