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

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

Doxygen comment added

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