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

Last change on this file since 2042 was 2042, checked in by suehring, 5 years ago

Bugfix, read and write restart data for wall fluxes; additional checks for wall fluxes

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