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

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

last commit documented

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