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

Last change on this file since 2084 was 2084, checked in by knoop, 5 years ago

Bugfix: Multigrid now usable with anelastic approximation

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