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

Last change on this file since 2024 was 2024, checked in by kanani, 6 years ago

changes related to urban surface model and output of ssws

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