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

Last change on this file since 1980 was 1973, checked in by maronga, 8 years ago

last commit documented

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