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

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

last commit documented

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