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

Last change on this file since 2312 was 2312, checked in by hoffmann, 4 years ago

various improvements of the LCM

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