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

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