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

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

Bugfix in definition of generic topography; missing check for microphysics added

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