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

Last change on this file since 2050 was 2050, checked in by gronemeier, 5 years ago

Implement turbulent outflow condition

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