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

Last change on this file since 2053 was 2053, checked in by gronemeier, 5 years ago

last commit documented

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