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

Last change on this file since 2026 was 2026, checked in by suehring, 8 years ago

Bugfix, enable output of s*2; calculation of domain_averaged perturbation energy; some formatting adjustments

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