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

Last change on this file since 2274 was 2274, checked in by Giersch, 7 years ago

Changed error messages

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