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

Last change on this file since 2263 was 2259, checked in by gronemeier, 7 years ago

Implemented synthetic turbulence generator

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