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

Last change on this file since 2011 was 2011, checked in by kanani, 5 years ago

changes related to steering and formating of urban surface model

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