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

Last change on this file since 2201 was 2201, checked in by suehring, 6 years ago

last commit documented

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