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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago

Forced header and separation lines into 80 columns

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