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

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago

last commit documented

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